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Foreword 


Statistics, mathematics and computer science have always been the favourite subjects in my 
academic career. The Ph.D. in automation and computer science engineering brought me to 
address challenging problems involving such disciplines. In particular, multi-object filtering 
concerns the joint detection and estimation of an unknown and possibly time-varying number 
of objects, along with their dynamic states, given a sequence of observation sets. Further, 
its distributed formulation also considers how to efficiently address such a problem over a 
heterogeneous sensor network in a fully distributed, scalable and computationally efficient 
way. Distributed multi-object filtering is strongly linked with statistics and mathematics 
for modeling and tackling the main issues in an elegant and rigorous way, while computer 
science is fundamental for implementing and testing the resulting algorithms. This topic 
poses significant challenges and is indeed an interesting area of research which has fascinated 
me during the whole Ph.D. period. 

This thesis is the result of the research work carried out at the University of Florence 
(Florence, Italy) during the years 2012-2014, of a scientific collaboration for the biennium 
2012-2013 with Selex ES (former SELEX SI, Rome, Italy) and of 6 months spent as a visiting 
Ph.D. scholar at the Curtin University of Technology (Perth, Australia) during the period 
January-July 2014. 
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Abstract 


The aim of the present dissertation is to address distributed tracking over a network of het¬ 
erogeneous and geographically dispersed nodes (or agents) with sensing, communication and 
processing capabilities. Tracking is carried out in the Bayesian framework and its extension 
to a distributed context is made possible via an information-theoretic approach to data fusion 
which exploits consensus algorithms and the notion of Kullback-Leibler Average (KLA) of 
the Probability Density Functions (PDFs) to be fused. 

The first step toward distributed tracking considers a single moving object. Consensus 
takes place in each agent for spreading information over the network so that each node can 
track the object. To achieve such a goal, consensus is carried out on the local single-object 
posterior distribution, which is the result of local data processing, in the Bayesian setting, 
exploiting the last available measurement about the object. Such an approach is called 
Consensus on Posteriors (CP). The first contribution of the present work [BCF14a] is an 
improvement to the CP algorithm, namely Parallel Consensus on Likelihoods and Priors 
(CLCP). The idea is to carry out, in parallel, a separate consensus for the novel information 
(likelihoods) and one for the prior information (priors). This parallel procedure is conceived 
to avoid underweighting the novel information during the fusion steps. The outcomes of the 
two consensuses are then combined to provide the fused posterior density. Furthermore, the 
case of a single highly-maneuvering object is addressed. To this end, the object is modeled 
as a jump Markovian system and the multiple model (MM) filtering approach is adopted for 
local estimation. Thus, the consensus algorithms needs to be re-designed to cope with this 
new scenario. The second contribution [BCF + 14b] has been to devise two novel consensus 
MM filters to be used for tracking a maneuvering object. The novel consensus-based MM 
filters are based on the First Order Generalized Pseudo-Bayesian (GPBi) and Interacting 
Multiple Model (IMM) filters. 

The next step is in the direction of distributed estimation of multiple moving objects. 
In order to model, in a rigorous and elegant way, a possibly time-varying number of objects 
present in a given area of interest, the Random Finite Set (RFS) formulation is adopted since 
it provides the notion of probability density for multi-object states that allows to directly ex¬ 
tend existing tools in distributed estimation to multi-object tracking. The multi-object Bayes 
filter proposed by Mahler is a theoretically grounded solution to recursive Bayesian tracking 
based on RFSs. However, the multi-object Bayes recursion, unlike the single-object coun- 
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terpart, is affected by combinatorial complexity and is, therefore, computationally infeasible 
except for very small-scale problems involving few objects and/or measurements. For this rea¬ 
son, the computationally tractable Probability Hypothesis Density (PHD) and Cardinalized 
PHD (CPHD) filtering approaches will be used as a first endeavour to distributed multi¬ 
object filtering. The third contribution [BCF + 13a] is the generalisation of the single-object 
KLA to the RFS framework, which is the theoretical fundamental step for developing a novel 
consensus algorithm based on CPHD filtering, namely the Consensus CPHD (CCPHD). Each 
tracking agent locally updates multi-object CPHD, i.e. the cardinality distribution and the 
PHD, exploiting the multi-object dynamics and the available local measurements, exchanges 
such information with communicating agents and then carries out a fusion step to combine 
the information from all neighboring agents. 

The last theoretical step of the present dissertation is toward distributed filtering with 
the further requirement of unique object identities. To this end the labeled RFS framework 
is adopted as it provides a tractable approach to the multi-object Bayesian recursion. The 5- 
GLMB filter is an exact closed-form solution to the multi-object Bayes recursion which jointly 
yields state and label (or trajectory) estimates in the presence of clutter, misdetections and as¬ 
sociation uncertainty. Due to the presence of explicit data associations in the 5-GLMB filter, 
the number of components in the posterior grows without bound in time. The fourth con¬ 
tribution of this thesis is an efficient approximation of the <j-GLMB filter [FVPV15], namely 
Marginalized 5-GLMB (M5-GLMB), which preserves key summary statistics (i.e. both the 
PHD and cardinality distribution) of the full labeled posterior. This approximation also facil¬ 
itates efficient multi-sensor tracking with detection-based measurements. Simulation results 
are presented to verify the proposed approach. Finally, distributed labeled multi-object track¬ 
ing over sensor networks is taken into account. The last contribution [FVV + 15] is a further 
generalization of the KLA to the labeled RFS framework, which enables the development 
of two novel consensus tracking filters, namely the Consensus Marginalized 5-Generalized 
Labeled Multi-Bernoulli (CM-d'GLMB) and the Consensus Labeled Multi-Bernoulli (CLMB) 
tracking filters. The proposed algorithms provide a fully distributed, scalable and computa¬ 
tionally efficient solution for multi-object tracking. 

Simulation experiments on challenging single-object or multi-object tracking scenarios 
confirm the effectiveness of the proposed contributions. 


6 



Introduction 



Recent advances in wireless sensor technology have led to the development of large networks 
consisting of radio-interconnected nodes (or agents) with sensing, communication and pro¬ 
cessing capabilities. Such a net-centric technology enables the building of a more complete 
picture of the environment, by combining information from individual nodes (usually with 
limited observability) in a way that is scalable (w.r.t. the number of nodes), flexible and reli¬ 
able (i.e. robust to failures). Getting these benefits calls for architectures in which individual 
agents can operate without knowledge of the information flow in the network. Thus, taking 
into account the above-mentioned considerations, Multi-Object Tracking (MOT) in sensor 
networks requires redesigning the architecture and algorithms to address the following issues: 

• lack of a central fusion node; 

• scalable processing with respect to the network size; 

• each node operates without knowledge of the network topology; 

• each node operates without knowledge of the dependence between its own information 
and the information received from other nodes. 

To combine limited information (usually due to low observability) from individual nodes, 
a suitable information fusion procedure is required to reconstruct, from the node information, 
the state of the objects present in the surrounding environment. The scalability requirement, 
the lack of a fusion center and knowledge on the network topology call for the adoption 
of a consensus approach to achieve a collective fusion over the network by iterating local 
fusion steps among neighboring nodes [OSFM07, XBL05, CA09, BC14]. In addition, due to 
the possible data incest problem in the presence of network loops that can causes double 
counting of information, robust (but suboptimal) fusion rules, such as the Chernoff fusion 
rule [CT12,CCM10] (that includes Covariance Intersection (Cl) [JU97, Jul08] and its gener¬ 
alization [MahOO]) are required. 

The focus of the present dissertation is distributed estimation , from the single object to 
the more challenging multiple object case. 
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In the context of Distribute Single-Object Filtering (DSOF), standard or Extended or 
Unscented Kalman filters are adopted as local estimators, the consensus involves a single 
Gaussian component per node, characterized by either the estimate-covariance or the infor¬ 
mation pair. Whenever multiple models are adopted for better describing the motion of the 
object in the tracking scenario, multiple Gaussian components per node arise and consensus 
has to be extended to this multicomponent setting. Clearly the presence of different Gaussian 
components related to different motion models of the same object or to different objects imply 
different issues and corresponding solution approaches that will be separately addressed. In 
this single-object setting, the main contributions in the present work are: 

i. the development of a novel consensus algorithm, namely Parallel Consensus on Like¬ 
lihoods and Priors (CLCP), that carries out, in parallel, a separate consensus for the 
novel information (likelihoods) and one for the prior information (priors); 

ii. two novel consensus MM filters to be used for tracking a maneuvering object, namely 
Distributed First Order Generalized Pseudo-Bayesian (DGPBi) and Distributed Inter¬ 
acting Multiple Model (DIMM) filters. 

Furthermore, Distribute Multi-Object Filtering (DMOF) is taken into account. To model 
a possibly time-varying number of objects present in a given area of interest in the presence 
of detection uncertainty and clutter, the Random Finite Set (RFS) approach is adopted. 
The RFS formulation provides the useful concept of probability density for multi-object states 
that allows to directly extend existing tools in distributed estimation to multi-object track¬ 
ing. Such a concept is not available in the MHT and JPDA approaches [Rei79, FS85, FS86, 
BSF88, BSL95, BP99]. However, the multi-object Bayes recursion, unlike the single-object 
counterpart, is affected by combinatorial complexity and is, therefore, computationally infea¬ 
sible except for very small-scale problems involving very few objects and/or measurements. 
For this reason, the computationally tractable Probability Hypothesis Density (PHD) and 
Cardinalized PHD (CPHD) filtering approaches will be used to address DMOF. It is recalled 
that the CPHD filter propagates in time the discrete distribution of the number of objects, 
called cardinality distribution, and the spatial distribution in the state space of such objects, 
represented by the PHD (or intensity function). It is worth to point out that there have been 
several interesting contributions [MahOO,CMC90,CJMR10,UJCR10,UCJ11] on multi-object 
fusion. More specifically, [CMC90] addressed the problem of optimal fusion in the case of 
known correlations while [MahOO, CJMR10, UJCR10, UCJ11] concentrated on robust fusion 
for the practically more relevant case of unknown correlations. In particular, [MahOO] first 
generalized Cl in the context of multi-object fusion. Subsequently, [CJMR10] specialized the 
Generalized Covariance Intersection (GCI) of [MahOO] to specific forms of the multi-object 
densities providing, in particular, GCI fusion of cardinality distributions and PHD functions. 
In [UJCR10], a Monte Carlo (particle) realization is proposed for the GCI fusion of PHD 
functions. The two key contributions in this thesis work are: 

i. the generalisation of the single-object KLA to the RFS framework; 



ii. a novel consensus CPHD (CCPHD) filter, based on a Gaussian Mixture (GM) imple¬ 
mentation. 

Multi-object tracking (MOT) involves the on-line estimation of an unknown and time- 
varying number of objects and their individual trajectories from sensor data [BP99, BSF88, 
Mah07b]. The key challenges in multi-object tracking also include data association uncer¬ 
tainty. Numerous multi-object tracking algorithms have been developed in the literature 
and most of these fall under the three major paradigms of: Multiple Hypotheses Tracking 
(MHT) [Rei79, BP99], Joint Probabilistic Data Association (JPDA) [BSF88] and Random 
Finite Set (RFS) filtering [Mah07b]. The proposed solutions are based on the recently in¬ 
troduced concept of labeled RFS that enables the estimation of multi-object trajectories in 
a principled manner [VV13]. In addition, labeled RFS-based trackers do not suffer from 
the so-called “ spooky effect ” [FSU09] that degrades performance in the presence of low de¬ 
tection probability like in the multi-object filters [VVC07, BCF + 13a, UCJ13]. Labeled RFS 
conjugate priors [VV13] have led to the development of a tractable analytic multi-object 
tracking solution called the 6-Generalized Labeled Multi-Bernoulli (d-GLMB) filter [VVP14]. 
The computational complexity of the <5-GLMB filter is mainly due to the presence of explicit 
data associations. For certain applications such as tracking with multiple sensors, partially 
observable measurements or decentralized estimation, the application of a 5-GLMB filter may 
not be possible due to limited computational resources. Thus, cheaper approximations to the 
d-GLMB filter are of practical significance in MOT. Core contribution of the present work is a 
new approximation of the 5-GLMB filter. The result is based on the approximation proposed 
in [PVV + 14] where it was shown that the more general Generalized Labeled Multi-Bernoulli 
(GLMB) distribution can be used to construct a principled approximation of an arbitrary 
labeled RFS density that matches the PHD and the cardinality distribution. The resulting 
filter is referred to as Marginalized S-GLMB (Md-GLMB) since it can be interpreted as a 
marginalization over the data associations. The proposed filter is, therefore, computationally 
cheaper than the 5-GLMB filter while preserving key summary statistics of the multi-object 
posterior. Importantly, the Md-GLMB filter facilitates tractable multi-sensor multi-object 
tracking. Unlike PHD/CPHD and multi-Bernoulli based filters, the proposed approximation 
accommodates statistical dependence between objects. An alternative derivation of the La¬ 
beled Multi-Bernoulli (LMB) filter [RVVD14] based on the newly proposed Md-GLMB filter 
is presented. 

Finally, Distributed MOT (DMOT) is taken into account. The proposed solutions are 
based on the above-mentioned labeled RFS framework that has led to the development of 
the d-GLMB tracking filter [VVP14]. However, it is not known if this filter is amenable to 
DMOT. Nonetheless, the M<5-GLMB and the LMB filters are two efficient approximations of 
the d-GLMB filter that 

• have an appealing mathematical formulation that facilitates an efficient and tractable 
closed-form fusion rule for DMOT; 
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• preserve key summary statistics of the full multi-object posterior. 

In this setting, the main contributions in the present work are: 

i. the development of the first distributed multi-object tracking algorithms based on the 
labeled RFS framework, generalizing the approach of [BCF + 13a] from moment-based 
filtering to tracking with labels; 

ii. the development of Consensus Marginalized 5-Generalized Labeled Multi-Bernoulli (CM5- 
GLMB) and Consensus Labeled Multi-Bernoulli (CLMB) tracking filters. 

Simulation experiments on challenging tracking scenarios confirm the effectiveness of the 
proposed contributions. 

The rest of the thesis is organized as follows. 

Chapter 2 - Background 

This chapter introduces notation, provides the necessary background on recursive Bayesian 
estimation, Random Finite Sets (RFSs), Bayesian multi-object filtering, distributed estima¬ 
tion and the network model. 

Chapter 3 - Distributed single-object filtering 

This chapter provides novel contributions on distributed nonlinear filtering with applications 
to nonlinear single-object tracking. In particular: i) a Parallel Consensus on Likelihoods and 
Priors (CLCP) filter is proposed to improve performance with respect to existing consensus 
approaches for distributed nonlinear estimation; n) a consensus-based multiple model filter 
for jump Markovian systems is presented and applied to tracking of a highly-maneuvering 
objcet. 

Chapter 4 - Distributed multi-object filtering 

This chapter introduces consensus multi-object information fusion according to an information- 
theoretic interpretation in terms of Kullback-Leibler averaging of multi-object distributions. 
Moreover, the Consensus Cardinalized Probability Hypothesis Density (CCPHD) filter is pre¬ 
sented and its performance is evaluated via simulation experiments. 

Chapter 5 - Centralized multi-object tracking 

In this chapter, two possible approximations of the 5-Generalized Labeled Multi-Bernoulli 
((5-GLMB) density are presented, namely i) the Marginalized 5-Generalized Labeled Multi- 
Bernoulli (M<5-GLMB) and n) the Labeled Multi-Bernoulli (LMB). Such densities will allow 
to develop a new centralized tracker and to establish a new theoretical connection to previous 
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work proposed in the literature. Performance of the new centralized tracker is evaluated via 
simulation experiments. 

Chapter 6 - Distributed multi-object tracking 

This chapter introduces the information fusion rules for Marginalized 5-Generalized Labeled 
Multi-Bernoulli (Md-GLMB) and Labeled Multi-Bernoulli (LMB) densities. An information- 
theoretic interpretation of such fusions, in terms of Kullback-Leibler averaging of labeled 
multi-object densities, is also established. Furthermore, the Consensus M5-GLMB (CM5- 
GLMB) and Consensus LMB (CLMB) tracking filters are presented as two new labeled dis¬ 
tributed multi-object trackers. Finally, the effectiveness of the proposed trackers is discussed 
via simulation experiments on realistic distributed multi-object tracking scenarios. 

Chapter 7 - Conclusions and future work 

The thesis ends with concluding remarks and perspectives for future work. 
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2.1 Network model 

Recent advances in wireless sensor technology has led to the development of large networks 
consisting of radio-interconnected nodes (or agents) with sensing, communication and pro¬ 
cessing capabilities. Such a net-centric technology enables the building of a more complete 
picture of the environment, by combining information from individual nodes (usually with 
limited observability) in a way that is scalable (w.r.t. the number of nodes), flexible and reli¬ 
able (i.e. robust to failures). Getting these benefits calls for architectures in which individual 
agents can operate without knowledge of the information flow in the network. Thus, taking 
into account the above-mentioned considerations, Object Tracking (OT) in sensor networks 
requires redesigning the architecture and algorithms to address the following issues: 
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• lack of a central fusion node; 

• scalable processing with respect to the network size; 

• each node operates without knowledge of the network topology; 

• each node operates without knowledge of the dependence between its own information 
and the information received from other nodes. 

The network considered in this work (depicted in Fig. 2.1) consists of two types of 
heterogeneous and geographically dispersed nodes (or agents): communication (COM) nodes 
have only processing and communication capabilities, i.e. they can process local data as well 
as exchange data with the neighboring nodes, while sensor (SEN) nodes have also sensing 
capabilities, i.e. they can sense data from the environment. Notice that, since COM nodes 
do not provide any additional information, their presence is needed only to improve network 
connectivity. 

Node Types 
SEN node 

COM node 

Node functions 

| Sensing 

Information 
exchange 

Local information 
Processing 

From a mathematical viewpoint, the network is described by a directed graph Q = (J\T, A) 
where Af = S U C is the set of nodes, S is the set of sensor and C the set of commu¬ 
nication nodes, and A C Af x AT is the set of arcs, representing links (or connections). 
In particular, (i,j) E A if node j can receive data from node i. For each node j G Af. 
A = {i G Af : (i,j) G .4,} denotes the set of in-neighbours (including j itself), i.e. the set of 
nodes from which node j can receive data. 

Each node performs local computation, exchanges data with the neighbors and gathers 
measurements of kinematic variables (e.g., angles, distances, Doppler shifts, etc.) relative to 
objects present in the surrounding environment (or surveillance area). The focus of this thesis 





Figure 2.1: Network model 
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will be the development of networked estimation algorithms that are scalable with respect 
to network size, and to allow each node to operate without knowledge of the dependence 
between its own information and the information from other nodes. 

2.2 Recursive Bayesian estimation 

The main interest of the present dissertation is estimation , which refers to inferring the val¬ 
ues of a set of unknown variables from information provided by a set of noisy measurements 
whose values depend on such unknown variables. Estimation theory dates back to the work 
of Gauss [Gau04] on determining the orbit of celestial bodies from their observations. These 
studies led to the technique known as Least Squares. Over centuries, many other techniques 
have been proposed in the field of estimation theory [Fisl2,Kol50,Str60,Tre04, Jaz07, AM12], 
e.g., the Maximum Likelihood , the Maximum a Posteriori and the Minimum Mean Square 
Error estimation. The Bayesian approach models the quantities to be estimated as ran¬ 
dom variables characterized by Probability Density Functions (PDFs), and provides an im¬ 
proved estimation of such quantities by conditioning the PDFs on the available noisy measure¬ 
ments. Hereinafter, we refer to the Bayesian approach as to recursive Bayesian estimation 
(or Bayesian filtering) , a renowned and well-established probabilistic approach for recursively 
propagating, in a principled way via a two-step procedure, a PDF of a given time-dependent 
variable of interest. The first key concept of the present work is, indeed, Bayesian filtering. 
The propagated PDF will be used to describe, in a probabilistic way, the behaviour of a 
moving object. In the following, a summary of the Bayes Filter (BF) is given, as well as 
a review of a well known closed-form solution of it, the Kalman Filter (KF) [Kal60,KB61] 
obtained in the linear Gaussian case. 

2.2.1 Notation 

The following notation is adopted throughout the thesis: col( where I is a finite set, de¬ 

notes the vector/matrix obtained by stacking the arguments on top of each other; diag( 
where I is a finite set, denotes the square diagonal matrix obtained by placing the arguments 
in the (i, i)-th position of the main diagonal; the standard inner product notation is denoted 
as 

(/, g) = J f{x) g(x)dx ; (2.1) 

vectors are represented by lowercase letters, e.g. x , x; spaces are represented by blackboard 
bold letters e.g. X, Y, L, etc. The superscript T stems for the transpose operator. 

2.2.2 Bayes filter 

Consider a discrete-time state-space representation for modelling a dynamical system. At 
each time k £ N, such a system is characterized by a state vector xj~ £ X C R nx , where n x is 
the dimension of the state vector. The state evolves according to the following discrete-time 
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stochastic model: 

Xk — i ( 2 . 2 ) 

where fk-i is a, possibly nonlinear, function; Wk-i is the process noise modeling uncertain¬ 
ties and disturbances in the object motion model. The time evolution (2.2) is equivalently 
represented by a Markov transition density 

Tk\k-i(x\C) , (2.3) 

which is the PDF associated to the transition from the state £ = Xk ~i to the new state 
x = x k . 

Likewise, at each time k, the dynamical system described with state vector Xk can be 

observed via a noisy measurement vector yk 6 Y C M. ny , where n y is the dimension of the 

observation vector. The measurement process can be modelled by the measurement equation 

Vk = hk(x k ,Vk) , (2.4) 


which provides an indirect observation of the state Xk affected by the measurement noise 
Vk- The modeling of the measurement vector is equivalently represented by the likelihood 
function 

9k(y\x ) , (2.5) 

which is the PDF associated to the generation of the measurement vector y = yk from the 
dynamical system with state x = Xk- 

The aim of recursive state estimation (or filtering) is to sequentially estimate over time 
Xk given the measurement history y\-k = {yi ,..., yk}- It is assumed that the PDF associated 
to yi-k given the state history x\-k = {aq,..., Xk] is 


k 

9l:k(yi:k\xi:k) = JJ 9 k{Vk\x k ) , 

K,= l 


( 2 . 6 ) 


i.e. the measurements y\-± are conditionally independent on the states aq ±. In the Bayesian 
framework, the entity of interest is the posterior density Pk{x) that contains all the informa¬ 
tion about the state vector xq given all the measurements up to time k. Such a PDF can 
be recursively propagated in time resorting to the well know Chapman-Kolmogorov equation 
and the Bayes’ rule [HL64] 


Pk\k-i{x) = J Vk\k-i(x\C) Pk-i(() d( , 


Pk(x) 


9k{Vk\x) p k \k-i(x) 

J 9k{Vk\C) Pk\k-i{()d( 


(2.7) 

( 2 . 8 ) 


given an initial density po{-)- The PDF Pk\k-i(’) is referred to as the predicted density, while 
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Pk(-) is the filtered density. 

Let us consider a multi-sensor centralized setting in which a sensor network (A7, A) conveys 
all the measurements to a central fusion node. Assuming that the measurements taken by 
the sensors are independent, the Bayesian filtering recursion can be naturally extended as 
follows: 


Pk\k-i{x) = j Tk\k-i(x\C) p k -i{0 d(, 

(2-9) 

II 9k(yl\x) Pk\k-i(x) 


, \ ieN 

Pk{x) = , . 

(2-10) 

/ II si (y l k\c) Pk\k-i(C)dC 

J A,- /v r 



ie N 


2.2.3 Kalman Filter 

The KF [Kal60, KB61, HL64] is a closed-form solution of (2.7)-(2.8) in the linear Gaussian 
case. That is, suppose that (2.2) and (2.4) are linear transformations of the state with 
additive Gaussian white noise, i.e. 


x k = A k _ iXk-x + w k _ 1 , (2.11) 

y k C k x k + v k , (2.12) 

where A k _\ is the n x x n x state transition matrix, C ' k is the n y x n x observation matrix, Wk-i 
and Vk are mutually independent zero-mean white Gaussian noises with covariances Qk-i 
and Rk, respectively. Thus, the Markov transition density and the likelihood functions are 

Tk\k-i(x\C) =M{x ; A fc _iC,<5fc-i) , (2.13) 

9k{y\x) =Af(y, C k x,R k ) , (2.14) 

where 

A7(x; TO, P) = |2t tP\~* e ~hU-m) T P~Hx-m) (2T5) 

is a Gaussian PDF. Finally, suppose that the prior density 

Pk-i(x) = Af(x; x fe _i,P fc _i) (2.16) 

is Gaussian with mean x k -i an d covariance P k -\ ■ Solving (2.7), the predicted density turns 
out to be 

Pk\k-i(x) =A f(x; x k \ k -i, P k \k-i) , (2-17) 

a Gaussian PDF with mean x k \k-i an d covariance P k \ k -i- Moreover, solving (2.8), the pos¬ 
terior density (or updated density ), turns out to be 

p k {x) =AT(x ; x k ,P k ) , (2.18) 
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i.e. a Gaussian PDF with mean x k and covariance P^. 

Remark 1. If the posterior distributions are in the same family as the prior probability 
distribution, the prior and posterior are called conjugate distributions, and the prior is called 
a conjugate prior for the likelihood function. The Gaussian distribution is a conjugate prior. 

The KF recursion for computing both predicted and updated pairs (x k -i, Pk-i) and 
(xk, Pk) is reported in Table 2.1. 

Table 2.1: The Kalman Filter (KF) 

for k = 1 , 2 ,... do 
Prediction 

Xk\k-i = -Afe-i^fe-i > Predicted mean 

Pk\k-i = + Qk- 1 > Predicted covariance matrix 


Correction 


Cfc — Vk CkXk |fc_i 

Sk = Rk + CkPk\k-iCj 

K k = P k \k-iCj Sj: 1 

Xk — x k\k— i T K k e k 

Pk = P k \ k -i ~ KkSkKj 

end for 

> Innovation 

> Innovation covariance matrix 
> Kalman gain 
o Updated mean 
> Updated covariance matrix 

In the centralized setting the network (A f, A) conveys 

all the measurements 

Vk = Ck x k + , 

(2.19) 

vi~Af(0,Ri) , 

(2.20) 

i £ J\f, to a fusion center in order to evaluate (2.10). The result amounts to stack all the 
information from all nodes i E J\T as follows 

Vk = col ( III ) 

(2.21) 


(2.22) 

R* = diag(^) 

(2.23) 


and then to perform the same steps of the KF. A summary of the Multi-Sensor KF (MSKF) 
is reported in Table 2.2. 

The KF has the advantage of being Bayesian optimal, but is not directly applicable to 
nonlinear state-space models. Two well known approximations have proven to be effective 
in situations where one or both the equations (2.2) and (2.4) are nonlinear: i) Extended KF 
(EKF) [MSS62] and n) Unscented KF (EKF) [JU97]. The EKF is a first order approximation 
of the Kalman filter based on local linearization. The UKF uses the sampling principles of 
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Table 2.2: The Multi-Sensor Kalman Filter (MSKF) 

for fc = 1,2,... do 

Prediction 

^k\k— 1 -A-k—lXk—1 

[> Predicted mean 

Pk\k-1 = M-iPk-iAl-i + Qk-i 

> Predicted covariance matrix 

Stacking 

y k = col(y|) ieAr 

C k = col (Cj) #6Ar 

R k = diag 

Correction 

&k Vk Ck%k\k— 1 

> Innovation 

Sk = Rk + CkPk\k-iCj 

i> Innovation covariance matrix 

Kk = Pk\k-iCjS ^ 

> Kalman gain 

%k %k\k— 1 P P-k^k 

> Updated mean 

Pk = Pk\k-1 - KkSkKj 

> Updated covariance matrix 

end for 



the Unscented Transform (UT) [JUDW95] to propagate the first and second order moments 
of the predicted and updated densities. 


2.2.4 The Extended Kalman Filter 

This subsection presents the EKF which is basically an extension of the linear KF whenever 
one or both the equations (2.2) and (2.4) are nonlinear transformations of the state with 
additive Gaussian white noise [MSS62], i.e. 

x k = fk- i(xfc_i) + w k - 1 , (2.24) 

Vk = hk(xk) + v k • (2.25) 


The prediction equations of the EKF are of the same form as the KF, with the transition 
matrix A^_\ of (2.11) evaluated via linearization about the updated mean Xk-i, he. 


Ak- i 


9fk- 1(') 


dx 


x=x k - 1 


(2.26) 


The correction equations of the EKF are also of the same form as the KF, with the observation 
matrix C of (2.12) evaluated via linearization about the predicted mean xuk- 1 , he. 


C k 


dhk(-) 


dx 


% &k\k — l 


(2.27) 


The EKF recursion for computing both predicted and updated pairs (xk- 1 , Pk-i) and ( Xk, Pk) 
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is reported in Table 2.3. 


Table 2.3: The Extended Kalman Filter (EKF) 


for k = 1,2,... do 


Prediction 


Ak- i — 


dfk- 1(- 


dx 


\x=x k -i 
•^k\k —1 Ak— l&k— 1 

i — Ak— \Pk— T Qk— 1 

Correction 








Uk 1 

Sk — -Rfc + 

ffc = ^Ifc-iC^ 1 
%k = 1 “ 1 “ K-k&k 

Pk = Pk\k-i ~ K k S k Kj 

end for 


> Linearization about the updated mean 

> Predicted mean 
> Predicted covariance matrix 


> Linearization about the predicted mean 

0 Innovation 
> Innovation covariance matrix 
o Kalman gain 
> Updated mean 
> Updated covariance matrix 


2.2.5 The Unscented Kalman Filter 

The UKF is based on the Unscented Transform (UT), a derivative-free technique capable 
of providing a more accurate statistical characterization of a random variable undergoing a 
nonlinear transformation [JU97]. In particular, the UT is a deterministic technique suited 
to provide an approximation of the mean and covariance matrix of a given random variable 
subjected to a nonlinear transformation via a minimal set of its samples. Let us consider 
the mean m and associated covariance matrix P of a generic random variable along with a 
nonlinear transformation function g(-), the UT proceeds as follows: 

• generates 2 n x + 1 samples X E ]^n a; x( 2 n :i; +i)^ ^ ie so ca ll ec t u-points, starting from the 
mean m with deviation given by the matrix square root E of P; 

• propagates the cr-points through the nonlinear transformation function </(•) resulting in 

Q g j^n^x^riaj+l). 


• calculates the new transformed mean rri' and associated covariance matrix P gg as well 
as the cross-covariance matrix Pxg of the initial and transformed cr-points. 

The pseudo-code of the UT is reported in Table 2.4. 

Given three parameters a a , f3 a and the weights c, w m and fU c are calculated exploiting 
the algorithm in Table 2.5. Moment matching properties and performance improvements are 
discussed in [JU97, WvdMOl] by resorting to specific values of a a , f3 a and n a . It is of common 
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Table 2.4: The Unscented Transformation (UT) 


procedure UT(m, P, g ) 

C, 1U C — UTW(o cr , [3 ( y 1 K(j) 

S = /p 


X = 


m. 


. m 


+ \fc. 0, E, —E 


G = g{X) 
m! = 


= GW C G T 
Pyg = GW c G T 
Return m'. P gg , P xg 

end procedure 


> Weights are calculated exploiting UTW in Table 2.5 

> 0 is a zero column vector 
[> g(-) is applied to each column of X 


Table 2.5: Unscented Transformation Weights (UTW) 


procedure UTW(a ff , f3 a , n a ) 

? = a 2 a (n x + K a ) - n x 

Wm = s(n x + ?) _1 

ivc 0) = ? (n x + ?) _1 + (1 - a 2 a + j3 a ) 

Wm'"’ 2nx \u!c 1, '"’ 2nx) = [ 2 ( 71 * + ?)]' 1 


Wm. = 
w c = 


w£\ .. .,w£ nx ) 


T 


T 


W c = (I - [w m ... w m \) diag( 

c = al(n x + K a ) 

Return c, w m , W c 
end procedure 


(o) 

Wc 


(I ~ [W m . ..W m }) 


practice to set these three parameters as constants, thus computing the weights once at the 
beginning of the estimation process. 

The UT can be applied in the KF recursion allowing to obtain a nonlinear recursive 
estimator known as UKF [JU97]. The pseudo-code of the UKF is shown in Table 2.6. 

The main advantages of the UKF approach are the following: 

• it does not require the calculation of the Jacobians (2.26) and (2.27). The UKF algo¬ 
rithm is, therefore, very suitable for highly nonlinear problems and represents a good 
trade-off between accuracy and numerical efficiency; 

• being derivative free, it can cope with functions with jumps and discontinuities; 

• it is capable of capturing higher order moments of nonlinear transformations [JU97]. 

Due to the above mentioned benefits, the UKF is herewith adopted as the nonlinear recursive 
estimator. 

The KF represents the basic tool for recursively estimating the state, in a Bayesian frame¬ 
work, of a moving object, i.e. to perform Single-Object Filtering (SOF) [FS85,FS86,BSF88, 
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Table 2.6: The Unscented Kalman Filter (UKF) 

for k = 1,2,... do 
Prediction 

•^k\k— 1> Pk\k —1 UT k—1 > —11 k—1 > /(')) 

-Pfc|fc-1 = ^fe|fe-l + Q 

Correction 

yk\k—li Cfc UT , P^|fc_i, hf)^j 

+ R 

%k &k\k—l T CkSfc (yVk — Vk\k-l) 

Pk = Pk\k-i - CkS^CZ 

end for 


BSL95,BSLK01,BP99]. A natural evolution of SOF is Multi-Object Tracking (MOT), which 
involves the on-line estimation of an unknown and (possibly) time-varying number of ob¬ 
jects and their individual trajectories from sensor data [BP99, BSF88, Mah07b]. The key 
challenges in multi-object tracking include detection uncertainty , clutter , and data associ¬ 
ation uncertainty. Numerous multi-object tracking algorithms have been developed in the 
literature and most of these fall under three major paradigms: Multiple Hypotheses Tracking 
(MHT) [Rei79, BP99]; Joint Probabilistic Data Association (JPDA) [BSF88]; and Random 
Finite Set (RFS) filtering [Mah07b]. In this thesis, the focus is on the RFS formulation since 
it provides the concept of probability density for multi-object state that allows to directly 
extend the single-object Bayesian recursion. Such a concept is not available in the MHT and 
JPDA approaches [Rei79,FS85,FS86,BSF88,BSL95,BP99]. 

2.3 Random finite set approach 

The second key concept of the present dissertation is the RFS [GMN97, Mah07b, Mah04, 
Mahl3] approach. In the case where the need is to recursively estimate the state of a possi¬ 
bly time varying number of multiple dynamical systems, RFSs allow to generalize standard 
Bayesian filtering to a unified framework. In particular, states and observations will be mod¬ 
elled as RFSs where not only the single state and observation are random, but also their 
number (set cardinality). The purpose of this section is to cover the aspects of the RFS 
approach that will be useful for the subsequent chapters. Overviews on RFSs and further 
advanced topics concerning point process theory, stochastic geometry and measure theory 
can be found in [Mat75, SKM95, GMN97, Mah03, VSD05, Mah07b]. 

Remark 2. It is worth pointing out that in multi-object scenarios there is a subtle difference 
between filtering and tracking. In particular, the first refers to estimating the state of a 
possibly time-varying number of objects without, however, uniquely identifying them, i.e. 
after having estimated a multi-object density a decision-making operation is needed to extract 
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the objects themselves. On the other hand, the term “tracking” refers to jointly estimating 
a possibly time-varying number of objects and to uniquely mark them over time so that 
no decision-making operation has to be carried out and object trajectories are well defined. 
Finally, it is clear that in single-object scenarios the terms filtering and tracking can be used 
interchangeably. 

2.3.1 Notation 

Throughout the thesis, finite sets are represented by uppercase letters, e.g. X, X. The 
following multi-object exponential notation is used 

h x = JJ h{x) , (2.28) 

xex 

where h is a real-valued function, with h 0 = 1 by convention [Mah07b]. The following 
generalized Kronecker delta [VV13,VVP14] is also adopted 

A f i, if x = y 

Sy(X) ± l , (2.29) 

I 0, otherwise 

along with the inclusion function, a generalization of the indicator function, defined as 

. f 1, if X C Y 

MX) 4 l . (2.30) 

I 0, otherwise 

The shortand notation 1 y(x) is used in place of ly({x}) whenever X = {x}. The cardinality 
(number of elements) of the finite set X is denoted by |X|. The following PDF notation will 
be also used. 

e~ x X n 

Poisson^] (n) = -y— , A € N , n G N, (2-31) 

! t—— , n€[a, b} 

b a ,oGM,6gM,o<6. (2.32) 

0, n(£[a,b ] 

2.3.2 Random finite sets 

In a typical multiple object scenario, the number of objects varies with time due to their 
appearance and disappearance. The sensor observations are affected by misdetection (e.g., 
occlusions, low radar cross section, etc.) and false alarms (e.g., observations from the envi¬ 
ronment, clutter, etc.). This is further compounded by association uncertainty, i.e. it is not 
known which object generated which measurement. The objective of multi-object filtering 
is to jointly estimate over time the number of objects and their states from the observation 
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history. 

In this thesis we adopt the RFS formulation, as it provides the concept of probability 
density of the multi-object state that allows us to directly generalize (single-object) estimation 
to the multi-object case. Indeed, from an estimation viewpoint, the multi-object system state 
is naturally represented as a finite set [WPS 10]. More concisely, suppose that at time k, 
there are Nf. objects with states Xkp, ..., a -k,N k , each taking values in a state space X C M n:c , 
i.e. the multi-object state at time k is the finite set 


N-k — {%k,h • • • ■ 'D'-.Y;. } C X. 


(2.33) 


Since the multi-object state is a finite set, the concept of RFS is required to model, in a 
probabilistic way, its uncertainty. 

An RFS A on a space X is a random variable taking values in A(X), the space of finite 
subsets of X. The notation J- n (JL) will be also used to refer to the space of finite subsets of 
X with exactly n elements. 

Definition 1. An RFS A is a random variable that takes values as (unordered) finite sets, 
i.e. a finite-set-valued random variable. 


At the fundamental level, like any other random variable, an RFS is described by its 
probability distribution or probability density. 


Remark 3. What distinguishes an RFS from a random vector is that: i) the number of 
points is random; n) the points themselves are random and unordered. 


The space A(X) does not inherit the usual Euclidean notion of integration and density. 
In this thesis, we use the Finite Set STatistics (FISST) notion of integration/density to 
characterize RFSs [Mah03,Mah07b]. 

From a probabilistic viewpoint, an RFS A is completely characterized by its multi-object 
density /(A). In fact, given /(A), the cardinality Probability Mass Function (PMF) p(n) 
that A have n > 0 elements and the joint conditional PDFs f(xi,X 2 , ■ ■ ■, x n \n) over X" given 
that A have n elements, can be obtained as follows: 


p{n) 

f(x i, • • • ,x n \n) 



.., x n }) dx 1 • • • dx n 


——- f({x l,...,Xn}) 
nl p(n) 


(2.34) 

(2.35) 


Remark 4. The multi-object density /(A) is nothing but the multi-object counterpart of 
the state PDF in the single-object case. 

In order to measure probability over subsets of X or compute expectations of random set 
variables, it is convenient to introduce the following definition of set integral for a generic 
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real-valued function g(X) (not necessarily a multi-object density) of an RFS variable X: 
r 00 1 r 

g(X)5X = y2— g({xi,...,x n })dxi---dx n (2.36) 

Jx frf n\ Jxn 

= g(0)+ [ g({x})dx+\[ g({xi, x 2 }) dxidx 2 H- 

Jx 2 Jx 2 

In particular, 


P(X) = Prob(X C X) = / f(X) SX 

Jx 


(2.37) 


measures the probability that the RFS X is included in the subset X of M. nx . The function 
P(X) is also known as the Belief-Mass Function (BMF) [Mah07b]. 

Remark 5. The BMF /3(X) is nothing but the multi-object counterpart of the state Cumu¬ 
lative Distribution Function (CDF) in the single-object case. 

It is also easy to see that, thanks to the set integral definition (2.36), the multi-object 
density, like the single-object state PDF, satisfies the trivial normalization constraint 


[ f(X) sx= 1. 
Jx 


(2.38) 


It is worth pointing out that the multi-object density, while completely characterizing 
an RFS, involves a combinatorial complexity; hence simpler, though incomplete, character¬ 
izations are usually adopted in order to keep the Multi-Object Filtering (MOF) problem 
computationally tractable. In this respect, the first-order moment of the multi-object den¬ 
sity, better known as Probability Hypothesis Density (PHD) or intensity function, has been 
found to be a very successful characterization [Mah04, Mah03, Mah07b]. In order to define 
the PHD function, let us introduce the number of elements of the RFS X which is given by 


where 


N x = 4>x(x) dx , 

Jx 


$x(x) = <**(£) 

«6.Y 


(2.39) 


(2.40) 


We would like to define the PHD function d(x) of X over the state space X so that the 
expected number of elements of X in X is obtained by integrating d(j over X, i.e. 


E[IVx] = / d(x) dx . 

Jx 


Since 


E [N X ] = J N x f(X ) 5X 


4>x(x) dx 


f(X) 6X 


(2.41) 

(2.42) 

(2.43) 
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<t>x{x) f(X) SX 


dx , 


comparing (2.41) with (2.44), it turns out that 

d(x) ± E[<j> x (x)} = j tx(x) f(X) 5X 

Without loss of generality, the PHD function can be expressed in the form 


d(x) = n s(x) 


(2.44) 


(2.45) 


(2.46) 


where 

OO 

n = E[n] = E[ra(X)] = ^ np{n) (2.47) 

n =0 


is the expected number of objects and s(-) is a single-object PDF, called location density , 
such that 



(2.48) 


It is worth to highlight that, in general, the PHD function d(-) and the cardinality PMF 
p(-) do not completely characterize the multi-object distribution. However, for specific RFSs 
defined in the next section, the characterization is complete. 


2.3.3 Common classes of RFS 

A review of the common RFS densities is provided [Mah07b] hereafter. 

Poisson RFS 

A Poisson RFS X on X is uniquely characterized by its intensity function d(-). The Poisson 

RFSs have the unique property that the distribution of the cardinality of X is Poisson with 

mean 

D = f d{x) dx , (2.49) 

Jx 

and for a given cardinality the elements of X are i.i.d. with probability density 

six) = . (2.50) 

The probability density of X can be written as 

tt(X) = e~ D d{x) (2.51) 

xex 

The Poisson RFS is traditionally described as characterizing no spatial interaction or complete 
spatial randomness in the following sense. For any collection of disjoint subsets Fj G X, j € N, 
it can be shown that the count functions |X n Bi\ are independent random variables which 
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are Poisson distributed with mean 


D Ph = / d(x) dx , 

JBj 


(2.52) 


and for a given number of points occurring in Bi the individual points are i.i.d. according to 

d(x) l Bi (x) 


The procedure in Table 2.7 illustrates how to generate a sample from a Poisson RFS. 


(2.53) 


Table 2.7: Sampling a Poisson RFS 

X = 0 

Sample n Poisson pj 
for i = 1 ,,n do 
Sample Xi n >-/ s(-) 

X = XU{x} 

end for 


Independent identically distributed cluster RFS 

An i.i.d. cluster RFS X on X is uniquely characterized by its cardinality distribution p(-) 
and matching intensity function d(-). The cardinality distribution must satisfy 

OO « 

D = 'S'' np(n ) = / d(x)dx, (2-54) 

but can otherwise be arbitrary, and for a given cardinality the elements of X are i.i.d. with 
probability density 

s(x) = ® (2.55) 

The probability density of an i.i.d. cluster RFS can be written as 

7T(X) = \X\\p(\X\) n s(x) (2.56) 

x&X 

Note that an i.i.d. cluster RFS essentially captures the spatial randomness of the Poisson 
RFS without the restriction of a Poisson cardinality distribution. The procedure in Table 2.8 
illustrates how a sample from an i.i.d. RFS is generated. 

Bernoulli RFS 

A Bernoulli RFS X on X has probability q = 1 — r of being empty, and probability r of being 
a singleton whose only element is distributed according to a probability density p defined 
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Table 2.8: Sampling an i.i.d. RFS 

X = 0 

Sample n ~ p(-) 

for i = 1,.... ,n do 

Sample Xj 

rSJ ') 

X = X U {x} 

end for 


on X. The cardinality distribution of a Bernoulli RFS is thus a Bernoulli distribution with 
parameter r. A Bernoulli RFS is completely described by the parameter pair ( r,p (•)). 

Multi-Bernoulli RFS 


A multi-Bernoulli RFS X on X is a union of a fixed number of independent Bernoulli RFSs 
with existence probability r® G (0,1) and probability density pW(-) defined on X for 
i = 1,..., I, i.e. 


A = (J A» 


(2.57) 


i =1 

It follows that the mean cardinality of a multi-Bernoulli RFS is 

i 






(2.58) 


i— 1 


A multi-Bernoulli RFS is thus completely described by the corresponding multi-Bernoulli 
parameter set j(r0),p0)(.)j j Its probability density is 


„ w =n(i-r») X 

3 =1 


(2.59) 


For convenience, probability densities of the form (2.59) are abbreviated by the form 7r(X) = 
{(r(*),p(*) j | A multi-Bernoulli RFS jointly characterizes unions of non-interacting points 

with less than unity probability of occurrence and arbitrary spatial distributions. The pro¬ 
cedure in Table 2.9 illustrates how a sample from a multi-Bernoulli RFS is generated. 


2.3.4 Bayesian multi-object filtering 

Let us now introduce the basic ingredients of the MOF problem [Mah07b, Mah04, Mahl3, 
Mahl4], i.e. the object RFS X^ c X at time k and the RFS Yi of measurements gathered 
by node i G M at time k. It is also convenient to define the overall measurement RFSs at 
time k, 

Yk = U (2.60) 

ieN 
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Table 2.9: Sampling a multi-Bernoulli RFS 

X = 0 

for * = 1,..., I do 

Sample u Uniforrriroi] 
if u < then 

Sample x ~ pW(-) 

X = XO{x} 
end if 
end for 


and up to time k, 


k 


Y uk = U • 

K,= l 


(2.61) 


In the random set framework, the Bayesian approach to MOF consists, therefore, of recur¬ 
sively estimating the object set Xy. conditioned to the observations Yy-k- The object set is 
assumed to evolve according to a multi-object dynamics 


Xk = &k-i(Xk-i) U Bk- 1 


(2.62) 


where Bk -1 is the RFS of new-born objects at time k — 1 and 

<f> fc -i(X) = (J (j) k - i(x) , (2.63) 

xex 

! {£c / } , with survival probability Psk-i 

(2.64) 

0, otherwise 

Notice that according to (2.62)-(2.64) each object in the set Xy. is either a new-born object 
from the set Bk -1 or an object survived from Xk- 1 , with probability Ps,k- 1 , and whose 
state vector x' has evolved according to the single-object dynamics (2.2) with x = Xk~\ and 
x' = Xk■ In a similar way, observations are assumed to be generated, at each node i G Af, 
according to the measurement model 


Yi = ¥ k {X k )OCi (2.65) 

where C l k is the clutter RFS (i.e. the set of measurements not due to objects) at time k and 
node i, and 


%(x) = u r k ( X ) , 

xex 


fc(*) 


{yk} , with detection probability Po.k 
0 , otherwise 


( 2 . 66 ) 

(2.67) 
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Notice that according to (2.65)-(2.67) each measurement in the set Y k is either a false one 
from the clutter set C l k or is related to an object in X k , with probability Pn,k-> according to 
the single-sensor measurement equation (2.4). Hence it is clear how one of the major benefits 
of the random set approach is to directly include in the MOF problem formulation (2.62)- 

(2.67) fundamental practical issues such as object birth and death, false alarms (clutter) 
and missed detections. Following a Bayesian approach, the aim is to propagate in time the 
posterior ( filtered ) multi-object densities Tr k {X) of the object set Xk given Y\. k , as well as the 
prior (predicted ) ones ^k\k-i(X) of X k given Y\ :k _i. Exploiting random set theory, it has 
been found [Mah07b] that the multi-object densities follow recursions that are conceptually 
analogous to the well-known ones of Bayesian nonlinear filtering, i.e. 

T k\k-l(X ) = 5Z , 

gk{Y k \z) Tr k \ k _ 1 {z)dz 

where (p k \k-i{ • | ■) i s the multi-object transition density to time k , g k ( ■ | •) is the multi-object 
likelihood function at time k [Mah03, Mah07a, Mah07b]. 

Remark 6. Despite this conceptual resemblance, however, the multi-object Bayes recursions 

(2.68) -(2.69), unlike their single-object counterparts, are affected by combinatorial complexity 
and are, therefore, computationally infeasible except for very small-scale MOF problems 
involving few objects and/or measurements. 



( 2 . 68 ) 

(2.69) 


For this reason, the computationally tractable PHD and Cardinalized PHD (CPHD) fil¬ 
tering approaches will be followed and briefly reviewed. 


2.3.5 CPHD filtering 

The CPHD filter [Mah07a, VVC07] propagates in time the cardinality PMFs p k \k-i( n ) and 
Pfc(n) as well as the PHD functions d k \ k _i(x) and d k (x) of X k given Y\- k -i and, respectively, 
Y\- k assuming that the clutter RFS, the predicted and filtered RFSs are i.i.d. cluster processes 
(see (2.56)). The resulting CPHD recursions (prediction and correction) are as follows 


Prediction 

n 

Pk\k-l{ n ) = X] Pb( n ~j) PS,k\k—lU) 

j=o (2.70) 

d k\k-ii x ) = db,k{x) + J (p k \k-ii x \0 Ps,k-i{() 4-i(C) d( 

Correction 
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Gk(d k \k~i(-) ,Y k ,n) p k \ k -i(n) 

Pk(n) = --- 

£ £°(4|fc-i(0 ,Y k ,'n) Pk\k-i(v) (2.71) 

r/=0 

dk(x ) = Gv k {x) 4| fe _i(x) 

where: Pb,k(') is the assumed birth cardinality PMF; Ps,k\k-i{') is the cardinality PMF of 
survived objects, given by 

00 ( t. \ 

PsMk-iU) = £1.1 P SM (1 - Ps,k) h ~ j Pk-i(t) ; (2.72) 

db,k(') is the assumed PHD function of new-born objects; Pk\k-i( x \C) is the state transition 
PDF (2.3) associated to the single object dynamics; the generalized likelihood functions 
Gki'i ■)') and Gy k (-) have cumbersome expressions which can be found in [VVC07]. 

Enforcing Poisson cardinality distributions for the clutter and object sets, the update 
of the cardinality PMF is no longer needed and the CPHD filter reduces to the PHD filter 
[Mah03, VM06]. The CPHD filter provides enhanced robustness with respect to misdetections, 
clutter as well as various uncertainties on the multi-object model (i.e. detection and survival 
probabilities, clutter and birth distributions) at the price of an increased computational load: 
a CPHD recursion has 0(m 3 n max ) complexity compared to 0(mn max ) of PHD, m being the 
number of measurements and n max the maximum number of objects. There are essentially 
two ways of implementing PHD/CPHD filters, namely the Particle Filter (PF) or Sequential 
Monte Carlo (SMC) [VSD05] and the Gaussian Mixture (GM) [VM06,VVC07] approaches, 
the latter being cheaper in terms of computation and memory requirements and, hence, by 
far preferable for distributed implementation on a sensor network wherein nodes have limited 
processing and communication capabilities. 

2.3.6 Labeled RFSs 

The aim of MOT is not only to recursively estimate the state of the objects and their time 
varying number, but also to keep track of their trajectories. To this end, the notion of label 
is introduced in the RFS approach [VV13,VVP14] so that each estimated object can be 
uniquely identified and its track reconstructed. 

Let C : XxL —> L be the projection C((x,£)) = £, where X C M na: , L = {cq : i G N} 
and ccj’s are distinct. To estimate the trajectories of the objects they need to be uniquely 
identified by a (unobserved) label drawn from a discrete countable space L = {ccj : i G N}, 
where the ccj’s are distinct. To incorporate object identity, a label £ e L is appended to 
the state x of each object and the multi-object state is regarded as a finite set on X x L, i.e. 
x = (x, £) 6 X x L. However, this idea alone is not enough since L is discrete and it is possible 
(with non-zero probability) that multiple objects have the same identity, i.e. be marked with 
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the same label. This problem can be alleviated using a special case of RFS called labeled 
RFS [VV13, VVP14], which is, in essence, a marked RFS with distinct labels. Then, a finite 
subset X^ of X x L has distinct labels if and only if X and its labels T{ -Xj = j T(x) : x £ X} 
have the same cardinality, i.e. |X| = |£(X)|. The function A(X) A <5|x|(|£(X)|) is called the 
distinct label indicator and assumes value 1 if and only if the condition |X| = |£(X)| holds. 

Definition 2. A labeled RFS with state space X and (discrete) label space L is an RFS on 
XxL such that each realization X has distinct labels, i.e. 

|X| = |£(X)| . (2.73) 

The unlabeled version of a labeled RFS with density 7r is distributed according to the 
marginal density [VV13] 

tt({.ti, ... ,x n }) = .,(x n ,£ n )}) . (2.74) 

The unlabeled version of a labeled RFS is its projection from XxL into X, and is obtained 
by simply discarding the labels. The cardinality distribution of a labeled RFS is the same as 
its unlabeled version [VV13]. 

Hereinafter, symbols for labeled states and their distributions are bolded to distinguish 
them from unlabeled ones, e.g. x, X, 7r, etc. 

2.3.7 Common classes of labeled RFS 

A review of the common labeled RFS densities is provided [VV13]. 

Labeled Poisson RFS 

A labeled Poisson RFS X with state space X and label space L = {aj : z £ N}, is a Poisson 
RFS A on X with intensity d(-), tagged with labels from L. A sample from such labeled 
Poisson RFS can be generated by the procedure reported in Table 2.10. 

Table 2.10: Sampling a labeled Poisson RFS 

X = 0 

Sample n Poisson pi 

for i = 1,... ,n do 
Sample x r*j d(-)/D 
X = X U {(x, cti)} 

end for 


The probability density of an LMB RFS is given by 

TT(X) = 4(|X|)(AX)) Poisson [D] (|X|) n ® , (2-75) 

xgX 
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where L (n) = {a* G ; 4( n )({4> •.. , £ n }) serves to check whether or not labels are dis¬ 
tinct. 

Remark 7. By a tracking point of view, each label of a labeled Poisson RFS X cannot 
directly refer to an object since all the corresponding states are sampled from the same 
intensity d(x), i.e. there is not a 1-to-l mapping between labels and location PDFs. Thus, it 
is not clear how to properly use such a labeled distribution in a MOT problem. 

Labeled independent identically distributed cluster RFS 

In the same fashion of the unlabeled RFSs, the labeled Poisson RFS can be generalized to 
the labeled i.i.d. cluster RFS by removing the Poisson assumption on the cardinality and 
specifying an arbitrary cardinality distribution. A sample from such a labeled i.i.d. RFS can 
be generated by the procedure reported in Table 2.11. 

Table 2.11: Sampling a labeled i.i.d. cluster RFS 

X = 0 

Sample n ~ p{ ) 
for i = 1 ,... ,n do 

Sample x r*j d(-)/D 
X = X U {(x, Oj)} 

end for _ 

The probability density of an LMB RFS is given by 

TT(X) = 4(1X1)(£(X)) p(|X|) n ® , (2-76) 

xex 

The same consideration drawn for the labeled Poisson RFS is inherited by the i.i.d. cluster 
RFS. 

Generalized Labeled Multi-Bernoulli RFS 

A Generalized Labeled Multi-Bernoulli (GLMB) RFS [VV13] is a labeled RFS with state 
space X and (discrete) label space L distributed according to 

7T(X) = A(X) 5>( c )(£(X)) [p®]* (2.77) 

ceC 

where: C is a discrete index set; w^ c \L) and j/ c ^ satisfy the normalization constraints: 

£5>( C >(L) = 1, (2.78) 

LCL cEC 

Jp( c \x ,£) dx = 1. (2.79) 
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A GLMB can be interpreted as a mixture of |C| multi-object exponentials 


ru (c) (£(X)) 


P 


0) 


1 x 


(2.80) 


Each term in the mixture (2.77) consists of the product of two factors: 

1. a weight uA c )(£(X)) that only depends on the labels £(X) of the multi-object state X; 


2. a multi-object exponential 


P 


(c) 


l x 


that depends on the entire multi-object state. 


The cardinality distribution of a GLMB is given by 


/>(»)= £ £>M(L), (2.81) 

LdTri L) CSC 

from which, in fact, summing over all possible n implies (2.78). The PHD of the unlabeled 
version of a GLMB is given by 


d(x) = EE^W) E 1 l(£)w^(L) . (2.82) 

ceCteL LCL 

Remark 8. The GLMB RFS distribution has the appealing feature of having a 1-to-l map¬ 
ping between labels and location PDFs provided by the single-object densities p( c \x,£), for 
each term of the multi-object exponential mixture indexed with c. However, in the context 
of MOT, it is not clear how to exploit (in particular how to implement) such a distribu¬ 
tion [VV13, VVP14], 


^-Generalized Labeled Multi-Bernoulli RFS 

A (^-Generalized Labeled Multi-Bernoulli (d-GLMB) RFS with state space X and (discrete) 
label space L is a special case of a GLMB with 


C = J’(L) x S, 

(2.83) 

w {c) (L) = w^(L) = w (I ^dj(L) , 

(2.84) 

P (c) (') =p (1 ’ a (-) = p^\-) , 

(2.85) 

i.e. is distributed according to 

7T(X) = A(X) £ wV#6i(£(X)) [p^l X 

(2.86) 

(/,06J(L)xS 

= A(X) £ i;(£(X))5;t«W)[p ( «] X , 

(2.87) 


ieH L) £e£ 


where 5 is a discrete space. In MOT a d-GLMB can be used (and implemented) to represent 
the multi-object densities over time. In particular: 
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• each finite set I £ <F(L) represents a different configurations of labels whose cardinality 
provides the number of objects; 

• each £ £ H represents a history of association maps, e.g. £ = (@i,... ,9k), where an 
association map at time 1 < j < k is a function which maps track labels at time j 
to a measurement at time j with the constraint that a track can generate at most one 
measurement, and a measurement can be assigned to at most one track; 

• the pair (/, £) is called hypothesis. 

To clarify the notion of 5-GLMB, let us consider the following two examples. 

El: Suppose the following two possibilities 

1. 0.4 chance of 1 object with label £\, i.e. I\ = {£i}, and density p(x,£ i); 

2. 0.6 chance of 2 objects with, respectively, labels £\ and £ 2 , i.e. I 2 = {^1,^2}, and 
densities p(x,£ 1) and p(x, £ 2 ). 

Then, the <j-GLMB representation is: 

tt(X) = 0.4ty l} (£(X)) P(x,£i) + 0.6fy liM (£(X)) p(x,h) p(x,£ 2 ) . (2.88) 


Notice that in this example there are no association histories, i.e. H = 0. Thus (2.88) 
has exactly two hypotheses (ii,0) and (I 2 ,0). 

E2: Let us now consider the previous example with S / 0: 

1. 0.4 chance of 1 object with label £\, i.e. I\ = {G}, association histories £1 and 
£2, i.e. there are two hypotheses (7i,£i) and (I\,£, 2 ), weights w^ 1 ’^ = 0.3 and 
w (h^ 2 ) — o.i ; densities p^\x,£ 1) and p^ 2 \x,£ 1); 

2. 0.6 chance of 2 objects with, respectively, label £\ and £ 2 , i.e. I 2 = {^1,^2}, associ¬ 
ation histories £1 and £2, i.e. there are two hypotheses (I 2 , £1) and (I 2 , £2), weights 
w (h,C) = 0.4 and w^ 2 ’^ = 0.2, densities p^\x,£ 1), p^\x,£ 2), p^ 2 \x,£ 1) and 
p^ 2 \x,£ 2 ). 


Then, the (5-GLMB representation is: 


tt(X) = fy l} (£(X)) [0.3p^(x,G) + 0.1;p ( ^(x,4) 

fy 1)M (£(X)) |o.4p^(x,4) p^ 1 \x,£ 2 ) + 0.2p^\x,£ 1 ) p^\x,£ 2 ) 1(2.89) 


->(£2)/ 


Thus (2.89) has exactly four hypotheses (/i,£i), (J\, £2); (^2,Ci) and (J 2 , £2)- 

The weight u/b£) represents the probability of hypothesis (I, £) and p^ (x, £) is the probability 
density of the kinematic state of track £ for the association map history £. 
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Labeled Multi-Bernoulli RFS 

A labeled multi-Bernoulli (LMB) RFS X with state space X, label space L and (finite) 
parameter set {(r^\p^) : £ G L}, is a multi-Bernoulli RFS on X augmented with labels 
corresponding to the successful (non-empty) Bernoulli components. In particular, is the 
existence probability and is the PDF on the state space X of the Bernoulli component 
with unique label < £ L. The procedure in Table 2.12 illustrates how a sample 
from a labeled multi-Bernoulli RFS is generated. 

Table 2.12: Sampling a labeled multi-Bernoulli RFS 

X = 0 

for £ G L do 

Sample u Uniform r 0 u 
if u < then 

Sample x ~ p^(-) 

X = XU{(x,£)j 

end if 
end for 


The probability density of an LMB RFS is given by 

tt(X) = A(X) w(£(X))p x (2.90) 


where 


w(L) = n 1 l( 0 n i 1 - r(e) ) > (2.91) 

t&V\L 

p^\x) = p(x,£). (2.92) 

For convenience, the shorthand notation 7 r = | {r^\p^^j j will be adopted for the density 
of an LMB RFS. The LMB is also a special case of the GLMB [VV13, VVP14, RVVD14] 
having C with a single element (thus the superscript is simply avoided) and 

p( c )(x,f) =p(x,£) = p^\x) , (2.93) 

w {c \L) = w(L) . (2.94) 

2.3.8 Bayesian multi-object tracking 

The labeled RFS paradigm [VV13, VVP14], along with the mathematical tools provided by 
FISST [Mah07b], allows to formalize in a rigorous and elegant way the multi-object Bayesian 
recursion for MOT. 

For MOT, the object label is an ordered pair of integers £ = ( k,i ), where k is the time 
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of birth and i E N is a unique index to distinguish objects born at the same time. The 
label space for objects born at time k is L*, = {k} x N. An object born at time k has state 
x e X x Lfc. Hence, the label space for objects at time k (including those born prior to k), 
denoted as Lo : fc, is constructed recursively by Lo : fc = Lo : fc-i U L/,. (note that Lo : fc-i and L*, 
are disjoint). A multi-object state X at time k is a finite subset of X x Lo,*,. 

Suppose that, at time k, there are objects with states x^i,..., ^k,N k , each taking 
values in the (labeled) state space X x L 0: fc, and M*. measurements yk,\, ..., Uk,M k each taking 
values in an observation space Y. The multi-object state and multi-object observation , at time 
k, [Mah03, Mah07b] are, respectively, the finite sets 

Xfc = {x fci i,...,x fciJVfc } ,. (2.95) 

Yk {Uk, 1 ) • • • ) Vk,Mk } • ( 2 . 96 ) 


Let 7Tfc(-) denote the multi-object filtering density at time k, and 7r fc | fc _ 1 (-) the multi-object 
prediction density to time k (formally 7tfc(-) and ^k\k-i{') should be written respectively 
as 7Tfc(-|Yo,..., Yfc_i, Yf.), and 7T fe i fc _ 1 (-|lo,..., Yfe_i), but for simplicity the dependence on 
past measurements is omitted). Then, the multi-object. Bayes recursion propagates 7r*, in 
time [Mah03, Mah07b] according to the following update and prediction 

"fclfc—i(X) = J ^.| fc _i(X|Z)7T fc _ 1 (Z)5Z, 

9k(Y k \X) 7Tfc|fe_i(X) 

Sk(Yk\Z) 7T fc | fc _ 1 (Z) SZ 

where <Pk\k-i (' I ') is the labeled multi-object transition density to time k, gk( -1 -) is the multi¬ 
object likelihood function at time k, and the integral is a set integral defined, for any function 
/ : F{X x L) —» R, by 

r oo , .. 

//(X)^X=V— V / f{{(xi,h ),..., (x n ,£ n )}) dxi ■ ■ ■ dx n . (2.99) 

/ —' Til z —' hvn 

The multi-object posterior density captures all information on the number of objects, and 
their states [Mah07b]. The multi-object likelihood function encapsulates the underlying mod¬ 
els for detections and false alarms while the multi-object transition density embeds the un¬ 
derlying models of motion, birth and death [Mah07b]. 



(2.97) 

(2.98) 


Remark 9. Eqs. (2.97)-(2.98) represent the multi-object counterpart of (2.7)-(2.8). 


Let us now consider a multi-sensor centralized setting in which a sensor network (J\f, A) 
conveys all the measurement to a central fusion node. Assuming that the measurements 
taken by the sensors are independent, the multi-object Bayesian filtering recursion can be 
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naturally extended as follows: 


^k\k— 1 (-^) J 

^ fe , fe _i(X|Z) 

7Tfc_i(Z) SZ , 

(2.100) 

7T *(X) = - 

ru(ni x 

ie N 

r 

) 7r fc|A:-l( X ) 

(2.101) 


n <?i(Vfci z ) 7r feife- 1 (z) dz 

i&N 


Remark 10. Eqs. (2.100)-(2.101) represents the multi-object counterpart of (2.9)-(2.10) 
which has been made possible thanks to the concept of RFS densities. 


In the work [VV13], the full labeled multi-object Bayesian recursion is derived for the 
general class of GLMB densities, while in [VVP14] the analytical implementation of the 
multi-object Bayesian recursion for the d-GLMB density is provided. Moreover, in the work 
[RVVD14], a multi-object Bayesian recursion for the LMB density, which turns out not to be 
in a closed form, is discussed and presented. 

For convenience, hereinafter, explicit reference to the time index k for label sets will be 
omitted by denoting L_ = Lo : fc-i, B = L&, L = L_ U B. 


2.4 Distributed information fusion 


The third and fourth core concepts of this dissertation are the Kullback-Leibler Average 
(KLA) and consensus, which are two mathematical tools for distributing information over 
sensor networks. 

To combine limited information from individual nodes, a suitable information fusion pro¬ 
cedure is required to reconstruct, from the information of the various local nodes, the state of 
the objects present in the surrounding environment. The scalability requirement, the lack of a 
fusion center and knowledge on the network topology (see section 2.1) dictate the adoption of 
a consensus approach to achieve a collective fusion over the network by iterating local fusion 
steps among neighboring nodes [OSFM07,XBL05,CA09,BC14], In addition, due to the data 
incest problem in the presence of network loops that causes double counting of information, 
robust (but suboptimal) fusion rules, such as the Chernoff fusion rule [CT12,CCM10] (that 
includes Covariance Intersection [JU97, Jul08] and its generalization [MahOO]) are required. 

In this chapter, the KLA and consensus will be separately dealt with, respectively, as a 
robust suboptimal fusion rule and as a technique to spread information, in a scalable way, 
throughout sensor networks. 
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2.4.1 Notation 

Given PDFs p, q and a scalar a > 0, the information fusion © and weighting operators 
© [BCF + 13a, BC14, BCF14a] are defined as follows 

(pe s )W*?^ (2.102) 

{P, Q) 

, “ 0 P)W *M' <2 ' 103) 

It can be checked that the fusion and weighting operators satisfy the following properties: 

P.A (p © q) © h = p © (q © h) = p® q@h 

P.B p ® q = q @ p 

P.C (a (5) © p = a © (/? © p) 

p.d 1 Qp = p 

p.e a © (p © q) = (a Op) © (a © q) 

P.F (a + /3) Qp = (a Qp) © (/3 © q) 

for any PDFs p and q. positive scalars a and (5. 


2.4.2 Kullback-Leibler average of PDFs 


A key ingredient for networked estimation is the capability to fuse in a consistent way PDFs 
of the quantity to be estimated provided by different nodes. In this respect, a sensible 
information-theoretic definition of fusion among PDFs is the Kullback-Leibler Average (KLA) 
[BCMP11,BC14] relying on the Kullback-Leibler Divergence (KLD). Given PDFs {p l (-)} ie jg- 
and relative weights oj 1 > 0 such that = ^ their weighted KLA p(-) is defined as 


where 


p(-) = arginf Vw’La'l (p || p l ) (2.104) 

^ ieN 

D K l{p || P l ) = Jp(x) log dx (2.105) 


denotes the KLD between the PDFs p(-) and p l (•). In [BC14] it is shown that the weighted 
KLA in (2.104) coincides with the normalized weighted geometric mean (NWGM) of the 
PDFs, i.e. 


p{x) 



~ © (w ? © p*(z)) 

iSAC 


(2.106) 
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where the latter equality follows from the properties of the operators © and ©. Note that in 
the unweighted KLA = 1/ |A/"|, i.e. 



(2.107) 


If all PDFs are Gaussian, i.e. p l (-) = A /"(•; m*,P l ), p(-) in (2.106) turns out to be Gaussian 
[BC14], i.e. p(-) =A/■(• ; m, P^j . In particular, defining the information (inverse covariance) 
matrix 

n = p- 1 (2.108) 

and information vector 

q = P~ l m, (2.109) 

associated to the Gaussian mean m and covariance P, one has 


( 2 . 110 ) 

i&Sf 

q = £ Nq 1 . (2.111) 

i(zAf 

which corresponds to the well known Covariance Intersection (Cl) fusion rule [JU97]. Hence 
the KLA of Gaussian PDFs is a Gaussian PDF whose information pair (D, q) is obtained 
by the weighted arithmetic mean of the information pairs (D\q l ) of the averaged Gaussian 
PDFs. 


2.4.3 Consensus algorithms 


In recent years, consensus algorithms have emerged as a powerful tool for distributed com¬ 
putation over networks [OSFM07, XBL05] and have been widely used in distributed parame¬ 
ter/state estimation algorithms [OS07,KT08,CS10,CCSZ08,CA09,SSS09,BCMP11,FFTS10, 
LJ12, HSH+12, BCF+13a, BC14, OHD14, BCF14a, BCF+14b, FBC+12, BCF+13b, BCF+13c, 
BCF + 13d,BCF + 14c]. In its basic form, a consensus algorithm can be seen as a technique for 
distributed averaging over a network; each agent aims to compute the collective average of 
a given quantity by iterative regional averages, where the terms “ collective" and lL regional v 
mean “over all network nodes ” and, respectively, u over neighboring nodes only'". Consensus 
can be exploited to develop scalable and reliable distributed fusion techniques. 

To this end, let us briefly introduce a prototypal average consensus problem. Let node 
i E Al be provided with an estimate 9 l of a given quantity of interest 6. The objective is to 
develop an algorithm that computes in a distributed way, in each node, the average 


6 


1 

Ni 




( 2 . 112 ) 
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To this end, let 9q = 9\ then a simple consensus algorithm takes the following iterative form: 

§\ +l = Y, w iJ 0j , Vi G Af (2.113) 

jeAf * 

where the consensus weights must satisfy the conditions 

5> iJ ' = l, VieAA (2.114) 

> 0, Vi, j eAf. (2.115) 


Notice from (2.113)-(2.115) that at a given consensus step the estimate in any node is com¬ 
puted as a convex combination of the estimates of the neighbors at the previous consensus 
step. In other words, the iteration (2.113) is simply a regional average computed in node 
i, the objective of consensus being convergence of such regional averages to the collective 
average (2.112). Important convergence properties, depending on the consensus weights, can 
be found in [OSFM07,XBL05]. For instance, let us denote by II the consensus matrix whose 
generic (*, j)-element coincides with the consensus weight a(if j (f Af l then u 1 ’ 3 is taken 
as 0). Then, if the consensus matrix II is primitive and doubly stochastic 1 , the consensus 
algorithm (2.119) asymptotically yields the average (2.112) in that 


lirn 0 l = 9 . 

£—>•00 


Vi G Af. 


(2.116) 


A necessary condition for the matrix II to be primitive is that the graph Q associated with 
the sensor network be strongly connected [CA09]. Moreover, in the case of an undirected 
graph Q, a possible choice ensuring convergence to the collective average is given by the 
so-called Metropolis weights [XBL05, CA09]. 


u 1 ’ 3 = 


1 


1 + max{|A/"*|, |A/V|} 


ieAfJ e A7\i / j 


(0 


t i,i 


= i- E 


UJ 


h3 


(2.117) 

(2.118) 




2.4.4 Consensus on posteriors 

The idea is to exploit consensus to reach a collective agreement (over the entire network), 
by having each node iteratively updating and passing its local information to neighbouring 
nodes. Such repeated local operations provide a mechanism for propagating information 
throughout the whole network. The idea is to perform, at each time instant and in each 
node of the network i 6 Af, a local correction step to find the local posterior PDF followed 
by consensus on such posterior PDFs to determine the collective unweighted KLA of the 
posterior densities p\. 

X A non-negative square matrix II is doubly stochastic if all its rows and columns sum up to 1. Further, it 
is primitive if there exists an integer m such that all the elements of IT n are strictly positive. 
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Suppose that at time k, each agent i starts with the posterior p\ as the initial iterate p\. 0 , 
and computes the l -th consensus iterate by 

Pk,l= ® (v tJ ®Pk,i- 1 ) ( 2 - 119 ) 

jeAf i 

where u l ’ J > 0, satisfying YhjeN 1 ujl ' 3 = 1> are the consensus weights relating agent i to nodes 
j G M l . Then, using the properties of the operators © and ©, it can be shown that [BC14] 

Pk,l= ® (^©-Pi) (2.120) 

j£Af 

where u]’ J is the (i,j)-th entry of PP, with II the (square) consensus matrix with (i,j)-th 
entry u l,J l (it is understood that pi, is omitted from the fusion whenever lu'C = 0). 
More importantly, it was shown in [OSFM07, XBL05] that if the consensus matrix II is 
primitive, (i.e. non-negative, and there exists an integer m such that II m is positive) and 
doubly stochastic (all rows and columns sum to 1), then for any i , j £ A7, one has 

<2 - m) 

In other words, at time k, if the consensus matrix is primitive then the consensus iterate of 
each node in the network “tends” to the collective unweighted KLA (2.107) of the posterior 
densities [CA09,BC14], 

The Consensus on Posteriors (CP) approach to Distributed SOF (DSOF) is summarized 
by the algorithm of Table 2.13 to be carried out at each sampling interval k in each node 
i £ M. In the linear Gaussian case, the CP algorithm involves a Kalman filter for local 
prediction and correction steps and Cl for the consensus steps. In this case, it has been 
been proved [BCMP11, BC14] that, under suitable assumptions, the CP algorithm guaran¬ 
tees mean-squared bounded estimation error in all network nodes for any number L > 1 of 
consensus iterations. In the nonlinear and/or non Gaussian case, the PDFs in the various 
steps can be approximated as Gaussian and the corresponding means-covariances updated 
via, e.g., EKF [MSS62] or UKF [JU04], 

Remark 11. Whenever the number of consensus steps L tends to infinity, consensus on 
posteriors is unable to recover the solution of the Bayesian filtering problem. In fact, in the CP 
approach, the novel information undergoing consensus combined with the prior information, 
is unavoidably underweighted. 
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Table 2.13: Consensus on Posteriors (CP) 


procedure CP(Node i, Time k ) 

Prediction 


Pk\k-l( x k) = 

1 <Pk\k-i( x k\xk-i) Pk-i(x k -i) dxk -1 

Local Correction 

pU x k\yi) = | 

(ai(yi\-) ©Pfc| fc -i(0) (*) > 

L do 

0 L 1 ’ 3 

Consensus 

Pkfi( X k) =pU 

for 1 = 1,..., 

Pk,l( x k ) = 


je Af i 

end for 


P k ( x k) = P%,L( x k) 

end procedure 
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This chapter presents two applications of distributed information fusion for single-object 
filtering [BCF14a,FBC + 12, BCF + 14b]. The first application presents an improvement of the 
CP approach, described in subsection 2.4.4, and is applied to track a single non-maneuvering 
object. The second application takes into account the possibility of the object of being highly 
maneuvering. Thus, a multiple-model filtering approach is used to devise a consensus-based 
algorithm capable of tracking a single highly-maneuvering object. The effectiveness of the 
proposed algorithms is demonstrated via simulation experiments on realistic scenarios. 


3.1 Consensus-based distributed filtering 

The approach proposed in this section is based on the idea of carrying out, in parallel, a 
separate consensus for the novel information (likelihoods) and one for the prior information 
(priors). This parallel procedure is conceived as an improvement of the CP approach to 
avoid underweighting the novel information during the fusion steps. The outcomes of the two 
consensuses are then combined to provide the fused posterior density. 

3.1.1 Parallel consensus on likelihoods and priors 

The proposed parallel Consensus on Likelihoods and Priors (CLCP) approach to DSOF is 
summarized by the algorithm of Table 3.1 to be carried out at each sampling interval k in 
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each node i G A f. Notice that in the correction step, a suitable positive weight p\ > 1 is 


Table 3.1: Consensus on Likelihoods and Priors (CLCP) pseudo-code 


procedure CLCP(Node i, Time k) 

Prediction 

P l k\k-i( x ) = J Tk\k-i{x\() Pk-i(x) d( 

Consensus 

*u«) = ( 6 * 

[ l, ieC 

Pk\k-i,o( x ) =Pk\k-i( x ) 

for l = 1,... ,L do 

9k,i( x ) = ® 

je N i 

Pk\k-iA x )= ® (w tJ QXifc-i.j-iC*)) 

jeN i 

end for 
Correction 

Pk(x) = Pi\k-1,AA 0 (& 0 9k,L (*)) 

end procedure 


applied to the outcome of the consensus on likelihoods, in order to possibly counteract the 
underweighting of novel information. To elaborate more on this issue, observe that each local 
posterior PDF p\{x) resulting from application of the CLCP algorithm turns out to be equal 
to 

Pk(x) = ( ® (“I 3 0 Pi\k-A x )) ) 0 ( ® [pWl 0 9k (v l k\x)) ] • (3.1) 

VieN / \jeS ) 

As it can be seen, p\(x ) is composed of two parts: i) a weighted geometric mean of the priors 
|fc_i(* T ); n ) a weighted combination of the likelihoods (f k (y\ \x). The fact that in the first part 
the weights u’jf sum up to one ensures that no double counting of the common information 
contained in the priors can occur [Jul08,BJA12]. As for the second part, care must be taken 
in the choice of the scalar weights p\. For instance, in order to avoid overweighting some 
of the likelihoods, it is important that, for any pair i.j. one has p\oo 1 ^ < 1. In this way, 
each component of independent information is replaced with a conservative approximation, 
thus ensuring the conservativeness of the overall fusion rule [BJA12]. On the other hand, the 
product p\oj 1 ^ should not be too small in order to avoid excessive underweighting. Based 
on these considerations and recalling the consensus property (2.121), different strategies for 
choosing p\ can be devised. 
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A reasonable choice would amount to letting 

Pk = m in 

jeS 

whenever at least one of the weights oj 1 £ j is different from zero. Such a choice would assign to 
p k the closest value to the true number of agents |A7| of the network. This strategy is easily 
applicable when L = 1 but, unfortunately, for L > 1 requires that each node of the network 
can compute U L . In most settings, this is not possible and alternative solutions must be 
adopted. 



For example, when L> 1 and the matrix II is primitive and doubly stochastic, one has 
that u l £ 3 ~ 1/|A7| for any i,j. Then, in this case, one can let 


pi = |AT|. (3.3) 

Such a choice has the appealing feature of giving rise to a distributed algorithm converging 
to the centralized one as L tends to infinity. However, when only a moderate number of 
consensus steps is performed, it leads to an overweighting of some likelihood components. 
Furthermore, the number of nodes |AA| might be unknown, in particular for a time-varying 
network. Notice that when such a choice is adopted, due to the properties of the fusion 0 and 
weighting © operators, the CLCP algorithm can be shown to be mathematically equivalent 
(in the sense that they would generate the same p\(x)) to the optimal distributed protocol 
of [OCIO]. 

An alternative solution is to exploit consensus so as to compute, in a distributed way, a 
normalization factor to improve the filter performance while preserving consistency of each 
local filter. For example, an estimate of the fraction |<S|/|A/"| of sensor nodes in the network 
can be computed via the consensus algorithm 


4,i = £ W<J l = (3.4) 

jeAf 1 

with the initialization b k 0 = 1 if i £ S, and b l k 0 = 0 otherwise. Then, the choice 



1, i^l,L = 0 

— , otherwise 

b k,L 


(3.5) 


has the desirable property of ensuring that p\oj l k < 1 for any i.j. Another positive feature 
of (3.5) is that no a priori knowledge of the network is required. 
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3.1.2 Approximate CLCP 


It must be pointed out that unfortunately the CLCP recursion of Table 3.1 does not admit an 
exact analytical solution except for the linear Gaussian case. Henceforth it will be assumed 
that the noises are Gaussian, i.e. ~ A7(0, Qk) and 4~ A/"(0, R k ), but that the system, i.e. 
fk(') and/or h k (-), might be nonlinear. The main issue is how to deal with nonlinear sensors 
in the consensus on likelihoods as for the other tasks (i.e. consensus on priors, correction, 
prediction) it is well known how to handle nonlinearities exploiting, e.g., EKF [MSS62] or 
UKF [JU04] or particle filters [RAG04], Under the Gaussian assumption, the local likelihoods 
take the form 

9k( x ) = 9k(y % k I*) = Vk - R k) (3-6) 

Whenever sensor i is linear, i.e. h\ (x) = C l k x, 

g{{x) oc e-^ xTsn i x ~ 2xTSq ti (3.7) 


where 

m = (ci) T (Rir'ci, 

Svi = (ci) T (RlCvi, 

Then it is clear that multiplying, or exponentiating by suitable weights, the likelihoods g k {-) is 
equivalent to adding, or multiplying by such weights, the corresponding and Sq' k defined 
in (3.8). Then, in the case of linear Gaussian sensors, the consensus on likelihoods reduces 
to a consensus on the information pairs defined in (3.8). For nonlinear sensors, a sensible 
approach seems therefore to approximate the nonlinear measurement function h\(-) by a 
linear affine one, i.e. 

( x ) — C k (x — x\\ k _^j + h l [x k (3.9) 

and then replace the measurement equation (2.4) with y\ = C\xk + v k for a suitably defined 
pseudo-measurement 

vi = yl — y l k\k-\ + cix^ k _ l (3.10) 

The sensor linearization (3.9) can be carried out for example by following the EKF paradigm. 
As an alternative, exploiting the unscented transform and the unscented information filter 
[Lee08], from x^ k _ 1 and the relative covariance P^ k _i, one can get the d-points j and 
relative weights ujj for j = 0,1,..., 2n (n = dirn(x)), and thus compute 


yl\k- 


i = £ y l k\k-v y l k\k -1 = h l (*fcffc-i) 

3=0 
2 n 

= E (: y l k\k-i - yjfeffc-i) (%-i - 


T 


3=0 


s~ii _ pV x I l pi \ 

— r k\k-l yk\k-l) 


-1 


(3.11) 

(3.12) 

(3.13) 
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Summarizing the above derivations, the proposed consensus-based distributed nonlinear filter 
is detailed in Table 3.2. 

The reason for performing in parallel consensus on priors and likelihoods is the counter¬ 
action of underweighting of the novel information contained in the likelihoods by using the 
weighting factor p\. Clearly, such a choice requires transmitting, at each consensus step, 
approximately twice the number of floating-point data with respect to other algorithms like 
CP. Notice that, depending on the packet size, this need not necessarily increase the packet 
rate. 


Table 3.2: Analytical implementation of Consensus on Likelihoods and Priors (CLCP) 


procedure CLCP(Node i, Time t) 


Consensus 


If i € S: 

If ieC: 
O i 

< 4 |fc-l ,0 


«4, 0 =(cj) T (i4 r'ci 
<«i,o=(ci) T (4) _1 a 
= o 


H,o = 0 
! ( F L-i) 

(*v 1 ) 

for i = 1, ..., L do 


-l 

-i 


l k\k-l 


Likelihood 


m,, = E 1 

jeAf* 

H,i= 

j£Af i 


Prior 


^k\k-l,l ~ X] 1,1—1 

is 

jeN' 

end for 


Correction 


= %k-i, L + pI m\, L 
<A = Qk\k-1,L + Pk H,l 


Prediction 


H = K)-',ii=(Slir 1 qi 

from x\, Pi compute i\ +m , P :_, , via UKF 
end procedure 
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3.1.3 A tracking case-study 


To evaluate performance of the proposed approach, the networked object tracking scenario of 
fig. 3.1 has been simulated. The network consists of 50 COM nodes and 10 SEN measuring the 


xIO 4 



- Link 

— -Surveillance area 
□ COM node 
• SEN node 
- Object trajectory 


Figure 3.1: Network and object trajectory. 


Time Of Arrival (TOA), i.e. the object-to-sensor distance. A linear white-noise acceleration 
model [BSLK01, p. 270], with 4-dimensional state consisting of positions and velocities 
along the coordinate axes, is assumed for the object motion; with this choice of the state, 
sensors are clearly nonlinear. Standard deviation of the measurement noise has been set to 
100 [m], the sampling interval to 5[s], the duration of the simulation to 1500 [s], i.e. 300 
samples. In the above scenario, 1000 Monte Carlo trials with independent measurement 
noise realizations have been performed. The Position Root Mean Square Error (PRMSE) 
averaged over time, Monte Carlo trials and nodes is reported in Table 3.3 for different values 
of L (number of consensus steps) and different consensus-based distributed nonlinear filters. 
In all simulations, Metropolis weights [CA09,XBL05] have been adopted for consensus. For 
each filter and each L two values are reported: the top one refers to the PRMSE computed 
over whole the simulation horizon, whereas the bottom one refers to the PRMSE computed 
after the transient period. Three approaches are compared in the UKF variants: the CLCP, 
CP (Consensus on Posteriors) [BCMP11,BC14] and CL (Consensus on Likelihoods) wherein 
consensus is carried out on the novel information only 1 . For the CLCP, three different 

1 In order to have a fair comparison, the CL algorithm is implemented as in Table 3.2 with the difference 
that no consensus on the priors is performed. 
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choices for the weights p\. are considered. From Table 3.3, it can be seen that the weight 
choice p\ = \Af\ provides the best performance only when L is sufficiently large, while for 
few consensus steps the other choices are preferable. It is worth pointing out that the CL 
approach, which becomes optimal as L —>■ oo, guarantees a bounded estimation error only 
when the number L of consensus steps is sufficiently high (in Table 3.3 ” indicates a 

divergent PRMSE). 


Table 3.3: Performance comparison 



CP 

CL 

CLCP 

Choice of p l k 



(3.5) 

m 

min jesil/^L 3 ) 

L = 1 

1733 

- 

1909 

- 

1645 

513 

- 

407 

- 

332 

L = 2 

1299 

- 

1182 

1621 

821 

423 

- 

257 

545 

226 

CO 

II 

1041 

- 

491 

736 

480 

389 

- 

231 

129 

198 

■'tf 

II 

936 

4269 

391 

455 

355 

370 

3238 

219 

118 

182 

LO 

II 

868 

977 

341 

283 

299 

358 

134 

210 

111 

170 

L = 15 

681 

222 

230 

138 

173 

318 

94 

182 

85 

122 


3.2 Consensus-based distributed multiple-model filtering 

The section addresses distributed state estimation of jump Markovian systems and its appli¬ 
cation to tracking of a maneuvering object by means of a network of heterogeneous sen¬ 
sors and communication nodes. It is well known that a single-model Kalman-like filter 
is ineffective for tracking a highly maneuvering object and that Multiple Model (MM) fil¬ 
ters [AF70, BBS88, BSLK01, LJ05] are by far superior for this purpose. The contributions 
provide novel consensus MM filters to be used for tracking maneuvering objects with sensor 
networks. Two novel consensus-based MM filters are presented. Simulation experiments in 
a tracking case-study, involving a strongly maneuvering object and a sensor network charac¬ 
terized by weak connectivity, demonstrate the superiority of the proposed filters with respect 
to existing solutions. 
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3.2.1 Notation 


The following new notations are adopted throughout the next sections. V c and V& denote the 
sets of PDFs over the continuous state space M n and, respectively, of PMFs over the discrete 
state space 7 Z, i.e. 


V c 


V d 


jp(-) : M n —>■ M 



1 and p(x) > 0, Vx E 



= 1 and p? > 0, Vj E 1Z 

jeTZ 


(3-14) 

(3.15) 


3.2.2 Bayesian multiple-model filtering 

Let us now focus on state estimation for the jump Markovian system 

Xk = f(m k , x k -\) + w k -i (3.16) 

y k = h(m k , x k ) + v k (3.17) 

where m k E IZ = {1,2,..., r-} denotes the system mode or discrete state; w k is the pro¬ 
cess noise with PDF p w (m k , •); v k is the measurement noise, independent of w k , with PDF 
Pvijnki •)• It is assumed that the system can operate in r possible modes, each mode j E IZ be¬ 
ing characterized by a mode-matched model with state-transition function /(j, •) and process 
noise PDF p w (j, •), measurement function h(j,j and measurement noise PDF p v (j,-). Fur¬ 
ther, mode transitions are modelled by means of a homogeneous Markov chain with suitable 
transition probabilities 


Pjt = prob (m k = j\m k -i = t) , j, f E IZ (3.18) 

whose possible time-dependence is omitted for notational simplicity. It is well known that 
a jump Markovian system (3.16)-(3.18) can effectively model the motion of a maneuvering 
object [AF70,BSLK01, section 11.6]. For instance a Nearly-Constant Velocity (NCV) model 
can be used to describe straight-line object motion while Coordinated Turn (CT) models with 
different angular speeds can describe object maneuvers. In alternative, models with different 
process noise covariances (small for straight-line motion and larger for object maneuvers) 
can be used with the same kinematic transition function /(•). For multimodal systems, 
the classical single-model filtering approach is clearly inadequate. To achieve better state 
estimation performance, the MM filtering approach [BSLK01, section 11.6] can be adopted. 
MM filters can provide, in principle, the Bayes-optimal solution for the jump Markov system 
state estimation problem by running in parallel r mode-matched Bayesian filters (one for 
each mode) on the same input measurements. Each Bayesian filter provides, at each time 
k, the conditional PDFs pP k | fc _ 1 (-) and of the state vector x k given the observations 

y k — {yi, ■ ■ •, Vk} and the mode hypothesis m k = j ■ Further, for each hypothesis m k = j, the 
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conditional modal probability = prob (m k = j\y T ) is also updated. In summary, Bayesian 
MM filtering amounts to representing and propagating in time information of the state vector 
of the jump Markovian system in terms of the r mode-matched PDFs p/(-) ,pj .\ r (-) G V c , 
j G 1Z, along with the PMF y, k \ T = col(yi/,| r j ^ G Vd- Then, from such distributions the 
overall state conditional PDF is obtained by means of the PDF mixture 

Pfc|fc-i(0 = Ti\ k -iPl\k-i(-) = Tk,\k-iPk\k-i(-) , (3-19) 

j&n 

Pk(-) = J2 Pk\kPk(-) = Pk\kPk {0 , (3-20) 

j&n 

where p k | fc _i(-) = col^^O)) . gW and Pk (0 = col (pi-(')) ■ 

The Bayesian MM filter is summarized in Table 3.4. 


It is worth to highlight, however, that the filter in Table 3.4 is not practically irnple- 
mentable even in the simplest linear Gaussian case, i.e. when i) the functions f(j, •) and 
h(j. •) are linear for all j G IZ; and n) the PDFs £>/,(•), p w (j , •) and p v (j, •) are Gaussian for 
all j G TZ. In fact, both the mode fusion and the mixing steps (see Table 3.4) produce, just 
at the first time instant k = 1, a Gaussian mixture of r components. When time advances, 
the number of Gaussian components grows exponentially with time k as r k . To avoid this 
exponential growth, the two most commonly used MM filter algorithms, i.e. First Order Gen¬ 
eralized Pseudo-Bayesian (GPBi) and Interacting Multiple Model (IMM), adopt a different 
strategy to keep all mode-matched PDFs Gaussian throughout the recursions. In particular, 
the GPBi algorithm (see Table 3.5) performs at each time k, before the prediction step, 
a re-initialization of all mode-matched PDFs with a single Gaussian component having the 
same mean and covariance of the fused Gaussian mixture. Conversely, the IMM algorithm 
(see Table 3.6) carries out a mixing procedure to suitably re-initialize each mode-matched 
PDF with a different Gaussian component. In the linear Gaussian case, the correction and 
prediction steps of MM filters are carried out, in an exact way, by a bank of mode-matched 
Kalman filters. Further, the likelihoods gi needed to correct the modal probabilities are 
evaluated as 



(3.21) 


where el and S k are the innovation and relative covariance provided by the KF matched to 
mode j. Whenever the functions f(j, •) and/or h(j, •) are nonlinear, even correction and/or 
prediction destroy the Gaussian form of the mode-matched PDFs. To preserve Gaussianity, 
as it is common practice in nonlinear filtering, the posterior PDF can be approximated 
as Gaussian by making use either of linearization of the system model around the current 
estimate [MSS62], i.e. EKF, or of the unscented transform [JU04], i.e. UKF. 
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3.2.3 Networked multiple model estimation via consensus 

A fundamental issue for networked estimation is to consistently fuse PDFs of the continuous 
quantity, or PMFs of the discrete quantity, to be estimated coming from different nodes. 
A suitable information-theoretic fusion among continuous probability distributions has been 
introduced in subsection 2.4.2, which can be straightforwardly applied in this context. 

Let us now turn the attention to the case of discrete variables. Given PMFs //, v 6 Vd 
their KLD is defined as 

Dkl (t-i || v) = ^ lo § ( 3 - 22 ) 

According to (2.104), the weighted KLA of the PMFs ji 1 € V c , i £ Af, is defined as 

Ji = arg inf V w l D KL (/r || /p) . (3.23) 

^ Vd ftu v ’ 

The following result holds. 

Theorem 1 (KLA of PMFs). 

The weighted KLA in (3.23) is given by 


i -1 = coll 


—7 ieN 

K = - 


ll of 

ie N _ 

e n (©)' 


(3.24) 


ten ieN 


Analogously to the case of continuous probability distributions, fusion and weighting 
operators can be introduced for PMFs /i.v E V,i and u > 0, as follows: 


3 3 

H © v = r) = col(rj j ) ■ n with rf = Ji U = /P © iP, Vj G 7 Z 

2Jrv 

ten 

(u j Y 

u © r = r] = col(rj j ). n with w = loQ rP , Vj e TZ 

E 0) 

ten 


(3.25) 


Thanks to the properties p.a-p.f of the operators © and © (see subsection 2.4.1), the KLA 
PMF in (3.23) can be expressed as 


R — © © R 1 

ieN 


(3.26) 


The collective fusion of PMFs (3.26) can be computed, in a distributed and scalable fashion, 


via consensus iterations 


R\ ~ © ( a; * J © (4-i 

je AT i 


(3.27) 


initialized from r 1 0 = r 1 and with consensus weights satisfying > 0 and J2jeN l = 1- 
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3.2.4 Distributed multiple-model algorithms 


Let us now address the problem of interest, i.e. DSOF for the jump Markovian system 
(3.16)-(3.17) over a sensor network (A f, A) of the type modeled in section 2.1. In this setting, 
measurements in (3.17) are provided by sensor nodes, i.e. 


h(- 


lIk = co\( yl) , 

v / ies 

= coi(h‘(-, -)) ieS 

.. __i 


Vk 


-- col (4) 




(3.28) 

(3.29) 

(3.30) 


where the sensor measurement noise i G S, is characterized by the PDF p v i(-:, •)• The main 
difficulty, in the distributed context, arises from the possible lack of complete observability 
from an individual node (e.g. a communication node or a sensor node measuring only angle 
or range or Doppler-frequency shift of a moving object) which makes impossible for such a 
node to reliably estimate the system mode as well as the continuous state on the sole grounds 
of local information. To get around this problem, consensus on both the discrete state PMF 
and either the mode-matched or the fused continuous state PDFs can be exploited in order 
to spread information throughout the network and thus guarantee observability in each node, 
provided that the network is strongly connected (i.e., for any pair of nodes i and j there exists 
a path from i to j and viceversa) and the jump Markovian system is collectively observable 
(i.e. observable from the whole set S of sensors). 

Let us assume that at sampling time k, before processing the new measurements yj., each 
node % G M be provided with the prior mode-matched PDFs p^ k _ 1 (-) = col (pfcjfc-i (')) 

along with the mode PMF /4|fc_i = c °l (Mfcjfc-i) which are the outcomes of the previous 
local and consensus computations up to time k — 1. Then, sensor nodes i £ 5 correct 
both mode-matched PDFs and mode PMF with the current local measurement y\. while 
communication nodes leave them unchanged; the PDFs and PMF obtained in this way are 
called local posteriors. At this point, consensus is needed in each node i 6 M to exchange 
local posteriors with the neighbours and regionally average them over the subnetwork A/"* of 
in-neighbours. The more consensus iterations are performed, the faster will be convergence 
of the regional KLA to the collective KLA at the price of higher communication cost and, 
consequently, higher energy consumption and lower network lifetime. Two possible consensus 
strategies can be applied: 

Consensus on the Fused PDF - first carry out consensus on the mode PMFs, then fuse 
the local mode-matched posterior PDFs over 7 Z using the modal probabilities resulting 
from consensus, and finally carry out consensus on the fused PDF; 


Consensus on Mode-Matched PDFs - carry out in parallel consensus on the mode PMFs 
and on the mode-matched PDFs and then perform the mode fusion using the modal 
probabilities and mode-matched PDFs resulting from consensus. 
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Notice that the first approach is cheaper in terms of communication as the nodes need to 
exchange just a single fused PDF, instead of multiple (r = \IZ\) mode-matched PDFs. Further, 
it seems the most reasonable choice to be used in combination with the GPBi approach, as 
only the fused PDF is needed in the re-initialization step. Conversely, the latter approach is 
mandatory for the distributed IMM filter. Recall, in fact, that the IMM filter re-initializes 
each mode-matched PDF with a suitable mixture of such PDFs, totally disregarding the 
fused PDF. Hence, to spread information about the continuous state through the network, it 
is necessary to apply consensus to all mode-matched PDFs. 

Summing up, two novel distributed MM filters are proposed: the Distributed GPB\ 
(DGPBi) algorithm of Table 3.7 which adopts the Consensus on the Fused PDF approach 
to limit data communication costs, and the Distributed IMM (DIMM) algorithm of Table 
3.8 which, conversely, adopts the Consensus on Mode-Matched PDFs approach. Since both 
DGPBi and DIMM filters, like their centralized counterparts, propagate Gaussian PDFs 
completely characterized by either the estimate-covariance pair (x, P) or the information 
pair ( q = P^x^Pl = P~ 1 ), the consensus iteration is simply carried as a weighted arithmetic 
average of the information pairs associated to such PDFs, see (2.110)-(2.111). 

3.2.5 Connection with existing approach 

It is worth to point out the differences of the proposed DGPBi and DIMM algorithms of 
Tables 3.7 and, respectively, 3.8 with respect to the distributed IMM algorithm of reference 
[LJ12]. In the latter, each node carries out consensus on the newly acquired information and 
then updates the local prior with the outcome of such a consensus. Conversely, in the DGPBi 
and DIMM algorithms, each node first updates its local prior with the local new information, 
and then consensus on the resulting local posteriors is carried out. To be more specific, recall 
first that, following the Gaussian approximation paradigm and representing Gaussian PDFs 
with information pairs, the local correction for each node i and mode j, takes the form 

i + <^4 (3-31) 

4 = 4-i + S & ’ ( 3 * 32 ) 

where SlYN , Sq f f denote the innovation terms due to the new measurements y\. In particular, 
for a linear Gaussian sensor model characterized by measurement function h l (j,x) = C l ’ J x 
and measurement noise PDF p v i(j,-) = A7(-; 0, the innovation terms in (3.31)-(3.32) 
are given by 


= 

(c iJ ) T 


1 c iJ , 

(3.33) 

II 

/ ..a 


-1 . 


(c l ’A 

(r 1 ’A 

Vk- 

(3.34) 


Further, for a nonlinear sensor i, approximate innovation terms can be obtained 

making use either of the EKF or of the unscented transform as shown in [LJ12]. 
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Now, the following observations can be made. 

• The algorithm in [LJ12], which will be indicated below with the acronym DIMM-CL 
(DIMM with Consensus on Likelihoods ), applies consensus to such innovation terms 
scaled by the number of nodes, i.e. to |A7| 5Ll l k and |A7| <5(///,_ ] . Conversely, the 
DIMM algorithm of Table 3.8 first performs the correction (3.31)-(3.32) and then applies 
consensus to the local posteriors (o)f, Q 1 ^, and similarly does the GPBi algorithm of 
Table 3.7 with the fused information pairs (f I\,q l k ). 

• A further difference between DIMM and DIMM-CL concerns the update of the discrete 
state PMF: the former applies consensus to the posterior mode probabilities n k J , while 
the latter applies it to the mode likelihoods g k J . 

A discussion on the relative merits/demerits of the consensus approaches adopted for 
DIMM and, respectively, DIMM-CL is in order. The main positive feature of DIMM-CL is 
that it converges to the centralized IMM as L —>• oo. However, it requires a minimum number 
of consensus iterations per sampling interval in order to avoid divergence of the estimation 
error. In fact, by looking at the outcome of CL 


sn 


ij _ 

k,L ~ 


i£Af 


it can be seen that a certain number of consensus iterations is needed so that the local 
information has time to spread through the network and the system becomes observable or, 
at least detectable, from the subset of sensors S k = jj G S : / oj. 

The proposed approach does not suffer from such limitations thanks to the fact that 
the whole posterior PDFs are combined. Indeed, as shown in [BCMP11, BC14], the single¬ 
model consensus Kalman filter based on CP ensures stability, for any L > 1 under the 
assumptions of system observability from the whole network and strong network connectivity, 
in that it provides a mean-square bounded estimation error in each node of the network. 
Although the stability results in [BCMP11, BC14] are not proved for jump Markov and/or 
nonlinear systems, the superiority of DIMM over DIMM-CL for tracking a maneuvering object 
when only a limited number of consensus iterations is performed is confirmed by simulation 
experiments, as it will be seen in the next section. With this respect, since the number 
of data transmissions, which primarily affect energy consumption, is clearly proportional to 
L, in many situations it is preferable for energy efficiency to perform only a few consensus 
steps, possibly a single one. A further advantage of CP over CL is that the former does 
not require any prior knowledge of the network topology as well as on the number of nodes 
in order to work properly. As a by-product, the Cl approach can also better cope with 
time-varying networks where, due to node join/leave and link disconnections, the network 
graph is changing in time. On the other hand, the proposed approach does not approach the 
centralized filter as L —> oo since the adopted fusion rule follows a cautious strategy so as to 
achieve robustness with respect to data incest [BC14]. 
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3.2.6 A tracking case-study 


To assess performance of the proposed distributed multiple-model algorithms described in 
section 3.2.4, the 2D tracking scenario of Fig. 3.2 is considered. From Fig. 3.2 notice 
that the object is highly maneuvering with speeds ranging in [0,300] [m/s] and acceleration 
magnitudes in [0,1.5] g , g being the gravity acceleration. The object state is denoted by 
x = [x, x, y, y] T where (x, y) and (x, y) represent the object Cartesian position and, 
respectively, velocity components. The sampling interval is T s = 5[s] and the total object 
navigation time is 640[s], corresponding to 129 sampling intervals. Specifically, r = 5 different 
Coordinated-Turn (CT) models [BSLK01, FS85, FS86] are used in the MM algorithms. All 
models have the state dynamics 
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(3.35) 


(3.36) 


for five different constant angular speeds u £ {—1, —0.5,0, 0.5,1} [°/s]. Notice, in particular, 
that, taking the limit for oj —> 0, the model corresponding to uj = 0 is nothing but the well 
known Discrete White Noise Acceleration (DWNA) model [FS85, BSLK01], The standard 
deviation of the process noise is taken as a w = 0.1 [m/s 2 ] for the DWNA (a; = 0) model and 
a w = 0.5 [m/s 2 ] for the other models. Jump probabilities (3.18), for the Markov chain, are 
chosen as follows: 
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(3.37) 


All the local mode-matched filters exploit the UKF [JU04], All simulation results have 
been obtained by averaging over 300 Monte Carlo trials. No prior information on the ini¬ 
tial object position is assumed, i.e. in all filters the object state is initialized in the cen- 
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ter of the surveillance region with zero velocity and associated covariance matrix equal to 


diagjlO 8 ,10 4 ,10 8 ,10 4 }. 




Figure 3.2: Object trajectory considered in the simulation experiments 

A surveillance area of 50 x 50 [km 2 ] is considered, wherein 4 communication nodes, 4 
bearing-only sensors measuring the object Direction of Arrival (DOA) and 4 TOA are de¬ 
ployed as shown in Fig. 3.3. Notice that the presence of communication nodes improves 
energy efficiency in that it allows to reduce average distance among nodes and, hence, the 
required transmission power of each individual node, but also hampers information diffusion 
across the network. This type of network represents, therefore, a valid benchmark to test the 
effectiveness of consensus state estimators. The following measurement functions characterize 
the DOA and TOA sensors: 


! Z[(x — x l ) + j (y — y 1 )] , if i is a DOA sensor 

__ (3-38) 

y (x — x 1 ) 2 + (y — y 1 ) 2 , if i is a TOA sensor 


where (x l , y l ) represents the known position of sensor i in Cartesian coordinates. The 
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standard deviation of DOA and TOA measurement noises are taken respectively as cfdoa = 
1 [°] and otoa = 100 [m\. The number of consensus steps used in the simulations ranges from 
L = 1 to L = 5, reflecting a special attention towards energy efficiency. 


x 10 4 



Figure 3.3: Network with 4 DOA, 4 TOA sensors and 4 communication nodes. 

The following MM algorithms have been compared: centralized GPBi and IMM, dis¬ 
tributed GPBi and IMM (i.e. DGPBi and DIMM) as well as the DIMM-CL algorithm of 
reference [LJ12], The performance of the various algorithms is measured by the PRMSE. 
Notice that averaging is carried out over time and Monte Carlo trials for the centralized 
algorithms, while further averaging over network nodes is applied for distributed algorithms. 
Table 3.9 summarizes performance of the various algorithms using the network of Fig. 3.3. 
Notice that the PRMSEs of DIMM-CL are not included in Table 3.9 since for all the con¬ 
sidered values of L such an algorithm exhibited a divergent behaviour. These results are 
actually consistent with the performance evaluation in [LJ12] where, considering circular 
networks with N > 4 nodes, a minimum number of L = 10 is required to track the object. 
Moreover, considering the network in Fig. 3.3 without communication nodes, L = 80 con¬ 
sensus steps are required for non-divergent behaviour of the DIMM-CL in [LJ12]. Notice 
that such a number of consensus steps might still be too large for practical applications, 
which is also stated in [LJ12]. As it can be seen from table 3.9, DIMM performs significantly 


60 














3.2. CONSENSUS-BASED DISTRIBUTED MULTIPLE-MODEL FILTERING 


better than DGPBi and, due to the presence of communication nodes, the distributed case 
exhibits significantly worse performance compared to the centralized case. On the other hand 
it is evident that, even in presence of communication nodes (bringing no information about 
the object position) and of a highly maneuvering object, distributed MM algorithms still 
work satisfactorily by distributing the little available information throughout the network by 
means of consensus. In fact, despite these difficulties, the object can still be tracked with 
reasonable confidence. Fig. 3.4 shows, for L = 5, the differences between two nodes that 
are placed almost at opposite positions in the network. It is noticeable that the main diffi¬ 
culty in tracking by communication nodes arises when a turn is approaching, especially, in 
this trajectory, during the last turn. Clearly, one should use a higher number of consensus 
steps L to allow data to reach two distant points in the network and, hence, to have similar 
estimated trajectory and mode probabilities. 



Figure 3.4: Estimated object trajectories of node 1 (TOA, green dash-dotted line) and node 
2 (TOA, red dashed line) in the sensor network of fig. 3.3. The real object trajectory is the 
blue continuous line. 
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Table 3.4: Centralized Bayesian Multiple-Model (CBMM) filter pseudo-code 


procedure CBMM(node i, time k) 
Prediction 
for mode j e 1Z do 

Pi\k-i( x ) = j Pw{j,x- /(j,C)) Pfc_i(C) d ( 
ten 

end for 


Correction 


for mode j eH do 

j . . Pv (j, Vk~h (j , x)) p fe ifc_i(a;) 

Pk( x ) = 7 - 

/ Pu (j, Vk ~ h(j, 0) Pk\k-i{C) d( 


9 k =P 


Pk = 


(yk\y k \ 


m k = J 


I J 


9k Pk\k-l 


end for 


Z 9kPk\k-i 
ten 


Mode fusion 


Pk(x) = Y PkPU x ) 

jen 


Mixing 

Vj, t e U : 


Z Pi* Pk 


ien 


Mj £ll: jp k [x) = Y Pk 3 Pk( x ) 
ten 


Re-init ializat ion 

Vj ell: p{(x) = pj,{x) 

end procedure 
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Table 3.5: Centralized First-Order Generalized Pseudo-Bayesian (CGBPi) filter pseudo-code 


procedure CGBPi(node i, time k ) 


Prediction 


for mode j £ 7 Z do 

d pi 


from x J k _ 1 ,P£_ 1 compute via EKF or UKF 

/ 4 i*-i = Y.P30A 

end for 


ten 


Correction 


for mode j £ do 

from ^ k \ k _ v Pk \ k _:1 compute x J k . P 3 k , e 3 k , S{ via EKF or UKF 


9 k = 


det ( 2 , 1 *jSlj 


o~Wk) ( 5 fe) e l 


j 9kTk\k-i 
Tk = 


end for 


E t t 

9k 9k\k-l 


Mode fusion 


j 


*k =Y.Uk- 

jen 

p k = Z / 4 

jen 

Vj £ TZ : x{ = x k , P 3 k = Pk 

end procedure 


P l + ~ 4) (®k - *fc) 
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Table 3.6: Centralized Interacting Multiple Model (CIMM) filter pseudo-code 


procedure CIMM(node i, time k) 


algorithm of Table 3.5 



As in the GPBi algorithm of Table 3.5 


As in the GPBi algorithm of Table 3.5 


W e K : 

2^ Pp Tk 

for mode j E 1Z do 


—j 

x k = 

E 

,Ai At 

P'k x k 





p! — 

E 

t\j 

Pk 

n 


t&t 



end for 





Re-init ializat ion 

Vj €IZ: x{ = x{, P 3 k = P J k 

end procedure 
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Table 3.7: Distributed GPBi (DGPBi) filter pseudo-code 


procedure DGPBi(node i, time k) 


Prediction 

for mode j € TZ do 

given Xfc-i an( i ^k -1 compute x l ^ k _ 1 and P^ k _ 1 via EKF or UKF 

/4’fLi = Y.PjtTk-i 

end for 


Correction 


for mode j G TZ do 

if i £ S then 


from x^ k _ l ,P k |[,_i compute xj~ , P k ’ J , e fe ’ J via EKF or UKF 


k\k 


1,1 

9 k = 


det 

l ,3 _ J,3,,h3 


(2u^) 


5 


i.i 1,1 1,1 | M I 

7T,o — 9k Tk\k-i {^teK 9k Tk\k-i\ 

else if z € C then 

pi,3 — pi,3 ,i,3 _ n h3 

'k\k-V r k ~ r k\k-V 9k,0 ~ 9k\k-l 


n-1 


■^k 

end if 

end for 


Consensus on modal probabilities 


for mode j 6 TZ do 
for I = 1,..., L do 

k 1 © Aifcfi-i 


9k,i - ® [ a 
iSA/' i 

end for 


2,7 M 

% = 9 k ,L 

end for 


Mode fusion 


x 


El,l *1,' 

9k x k 

jeii 


n fi — E 9 k 

je tz 


pi,3 _1_ ( A* _ Ai,j\ (p. _ 

-Ufc ' y^kfl x k ) y^kfl x k J 


Consensus on the fused PDF 


(-ffc,o) > 4,0 “ ^k,0 x k,0 


O* — 

“fc ,0 — 

for J — I;,... /. do 

= E 

*eA/' i 

4,1 = E 1 

iSA/" 4 

end for 


Re-initialization 

(«li) 


Vj G TZ : xj? = 


-l 


end procedure 


«U. ff = (nil) 


-1 
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Table 3.8: Distributed IMM (DIMM) filter pseudocode 


procedure DIMM(node i, time k) 

I Prediction I 


As in the DGPBi algorithm of Table 3.7 


Correction 


for mode j £ IZ do 

if i € S then 


from x^ k _ v P^ k _ x compute x£ 0 ,P k % e'/, S k J via EKF or UKF 


1,1 

9k = 


det (2uS 1 ^ 


f e -H4 J )' (4 J ) M 


e 2 \~k j v“fc i ~k 
i.j i.j i.j i-X 1 — ^ 

Pk ,o — 9k Pk\k-i [^teTZ 9k Tk\k-i 

end if 

if i € C then 

pij _ pi,3 .i,3 _ ,i,3 

x k, 0 — X k\k-V r k,0 ~ r k\k-V Pk, 0 — Pk\k -1 

end if 

end for 

Parallel consensus on modal probabilities &: mode-matched PDFs 


for mode j € TZ do 

c,i,j _ (poY 1 n iJ - n iJ r iJ 
14 fc,0 — C k ’°J ’ — k,0 x k,0 

for l = 1,..., L do 

/4’fz = ® 0 Mfcfi-i 

leN* 

^k,i = Ct7 *’ 1 1» 9k,l = ul, ‘9kj -1 

*eA/' i iSA/" 4 

end for 

ij i.j 
Pk = Pk,L 

end for 

Mode fusion 


Vi€K:P‘ ; i=(n«)T4i = F«4i 

4 = E 4 J 4 i 

j£K 


n =e 

jen 


r k,L + { X k ~ x k,L) \ X k ~ x k. 


Mixing 

for mode j € 1Z do 

i,t\j _ 


VtSK : = PjtPk PjjPk 


,3,h 


-1 


~*,7 \ "" i,t 7 ~i ,7 

**r = Z. Pk x k, 

t&TL 

A 1 = E »P 

ten 

end for 

end procedure 


/ ~i,t ~i,t \ ( ~i,t. ~i,t \ 
[ x k ~ x k,L) [ x k - x k,L) 


P d.t . / ^L.L ^' 

fe ,L + ( 30 - * 


T 
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Table 3.9: Performance comparison for the sensor network of fig. 3.3 



GPBi 

IMM 

PRMSE [m\ 

538 

446 


PRMSE [m] 

DGPBi 

DIMM 

L = 1 

1603 

960 

L = 2 

1312 

781 

tr-i 

II 

CO 

183 

711 

II 

^4 

1103 

672 

lO 

II 

^4 

1052 

648 
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The focus of the next sections is on Distributed MOF (DMOF) over the network of section 
2.1 wherein each node (tracking agent) locally updates multi-object information exploiting the 
multi-object dynamics and the available local measurements, exchanges such information with 
communicating agents and then carries out a fusion step in order to combine the information 
from all neighboring agents. Specifically, a CPHD filtering approach (see subsection 2.3.5) 
to DMOF will be adopted [BCF + 13a]. Hence, the agents will locally update and fuse the 
cardinality PMF p(-) and the location PDF s(-) (see (2.56)) that, for the sake of brevity, will 
also be referred to as the CPHD. 

More formally, the DMOF problem over the network (AA, 4) can be stated as follows. 
Each node i 6 A must estimate at each time k £ {1,2,...} the CPHD of the unknown 
multi-object set X ^ in (2.62)-(2.67) given local measurements Yj? for all n > 1 up to time k 
and data received from all adjacent nodes j £ A F\ {?'} so that the estimated pair (/?}(•) , s},(-)) 
be as close as possible to the one that would be provided by a centralized CPHD filter 
simultaneously processing information from all nodes. 

Clearly local CPHD filtering is the same as the centralized one of (2.70)-(2.71), operating 
on the local measurement set Yj! instead of Yj,. Conversely, CPHD fusion deserves special 
attention and will be tackled in the next sections. 


4.1 Multi-object information fusion via consensus 

Suppose that, in each node i of the sensor network, a multi-object density P(X) is available 
which has been computed on the basis of the information (collected locally or propagated 
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from other nodes) which has become available to node %. Aim of this section is to investigate 
whether it is possible to devise a suitable distributed algorithm guaranteeing that all the 
nodes of the network reach an agreement regarding the multi-object density of the unknown 
multi-object set. An important feature the devised algorithm should enjoy is scalability, i.e. 
the processing load of each node must be independent of the network size. For this reason, 
approaches based on the optimal (Bayes) fusion rule [CMC90, MahOO] are ruled out. In 
fact, they would require, for each pair of adjacent nodes (i,j), the knowledge of the multi¬ 
object density f{X\I l n P) conditioned to the common information POP and, in a practical 
network, it is impossible to keep track of such a common information in a scalable way. Hence, 
some robust suboptimal fusion technique has to be devised. 


4.1.1 Multi-object Kullback-Leibler average 

The first important issue to be addressed is how to define the average of the local multi¬ 
object densities p(X). To this end, taking into account the benefit provided by the RFS 
approach about defining multi-object densities, it is possible to extend the notion of KLA of 
single-object PDFs introduced in subsection 2.4.2 to multi-object ones. Let us first define the 
notion of KLD to multi-object densities f(X ) and g(X) by 

DklU II g) = j f(X) !og(^^y) SX (4.1) 


where the integral in (4.1) must be interpreted as a set integral according to the definition 
(2.36). Then, the weighted KLA f kla(X ) of the multi-object densities p(X) is defined as 
follows 

fi<LA = arginf ^ u l D KL (f || f'j . (4.2) 

* i&N 


with weights c o l satisfying 


a/>0, £V = 1. 

ie M 


(4.3) 


Notice from (4.2) that the weighted KLA of the agent densities is the one that minimizes the 
weighted sum of distances from such densities. In particular, if N is the number of agents and 
a; 1 = 1/N for i = 1 ,N, (4.2) provides the (unweighted) KLA which averages the agent 
densities giving to all of them the same level of confidence. An interesting interpretation of 
such a notion can be given recalling that, in Bayesian statistics, the KLD (4.1) can be seen 
as the information gain achieved when moving from a prior g(X) to a posterior f(X). Thus, 
according to (4.2), the average PDF is the one that minimizes the sum of the information 
gains from the initial multi-object densities. Thus, this choice is coherent with the Principle 
of Minimum Discrimination Information (PMDI) according to which the probability density 
which best represents the current state of knowledge is the one which produces an information 
gain as small as possible (see [Cam70, Aka98] for a discussion on such a principle and its 
relation with Gauss’ principle and maximum likelihood estimation) or, in other words [Jay03]: 
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“the probability assignment which most honestly describes what we know should 
be the most conservative assignment in the sense that it does not permit one to 
draw any conclusions not warranted by the data”. 

The adherence to the PMDI is important in order to counteract the so-called data incest 
phenomenon, i.e. the unaware reuse of the same piece of information due to the presence of 
loops within the network. 

The following result holds. 

Theorem 2 (Multi-object KLA). 

The weighted KLA defined in (4-2) turns out to be given by 

n [rwf 

fKLA ( x ) = - -J - (4.4) 

J i&N 


4.1.2 Consensus-based multi-object filtering 

The identity (4.4) with u/ = 1/|A7| can be rewritten as 

fKLA(X) = ® . (4.5) 

Then, the global (collective) KLA (4.5), which would require all the local multi-object densi¬ 
ties to be available, can be computed in a distributed and scalabale way by iterating regional 
averages through the consensus algorithm 

fl+ i(X) = 0 © f!(xj) , Vi G M (4.6) 

je AT* 

with fo(X) = f l (X). In fact, thanks to properties P.A - p.f, it can be seen that 1 

fi(X) = © © f j W) - Vi G A I (4.7) 

j&N 

where u^ 3 is defined as the element (i, j) of the matrix f V. With this respect, recall that when 
the consensus weights are chosen so as to ensure that the matrix Ll is doubly stochastic, 
one has 

VTooT ©. < 4 -s) 


1 Notice that the scalar multiplication operator is defined only for strictly positive scalars. However, in 
equation (4.7) it is admitted that some of the scalar weights w\' ] be zero. The understanding is that, whenever 
uj is equal to zero, the corresponding multi-object density / J (K) is omitted from the addition. This can 
always be done, since for each i £ Af and each l, at least one of the weights uj \’ 3 is strictly positive. 
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Hence, as the number of consensus steps increases, each local multi-object density “tends” to 
the global KLA (4.5). 


4.1.3 Connection with existing approaches 

It is important to point out that the fusion rule (4.4), which has been derived as KLA of 
the local multi-object densities, coincides with the Chernoff fusion [CT12,CCM10] known as 
Generalized Covariance Intersection (GCI) for multi-object fusion, first proposed by Mahler 
[MahOO]. In particular, it is clear that / kla{X ) in (4.4) is nothing but the NWGM of 
the agent multi-object densities /®(X); it is also called Exponential Mixture Density (EMD) 
[UJCRIO.UCJll]. 

Remark 12. The name GCI stems from the fact that (4.4) is the multi-object counterpart 
of the analogous fusion rule for (single-object) PDFs [MahOO, Hur02, JBU06, Jul08] which, in 
turn, is a generalization of Cl originally conceived [JU97] for Gaussian PDFs. 


Remind that, given estimates X{ of the same quantity x from multiple estimators with relative 
covariances Pi and unknown correlations, their Cl fusion is given by 


P = 


-l 




-l 




X = PJ2 ^P^Xi 

ieN 


(4.9) 


The peculiarity of (4.9) is that, for any choice of the weights ad satisfying (4.3) and provided 
that all estimates are consistent in the sense that 


E 


(x - Xi) (x - Xi) 


< Pi 


Vi 


(4.10) 


then the fused estimate also turns out to be consistent, i.e. 


E 


(x — x) (x — x) 


< p. 


(4.11) 


It can easily be shown that, for normally distributed estimates, (4.9) is equivalent to 


p(x) = 



(4.12) 


where p i (■) = Af(-; Xi, Pi) is the Gaussian PDF with mean X{ and covariance Pi. This sug¬ 
gested to use (4.12) for arbitrary, possibly non Gaussian, PDFs. As a final remark, notice 
that the consistency property, which was the primary motivation that led to the development 
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of the Cl fusion rule, is in accordance with the PMDI discussed in Section 4.1.1, which indeed 
represents one of the main positive features of the considered multi-object fusion rule. 


4.1.4 CPHD fusion 

Whenever the object set is modelled as an i.i.d. cluster process, the agent multi-object 
densities to be fused take the form 


f l (X) = \X\\ p\\X\) n^(x) (4.13) 

x£X 


where ( p l {n ) ,s 1 (x)) is the CPHD of agent i. In [CJMR10] it is shown that in this case the 
GCI fusion (4.4) yields 

j(x)=\x\\p(\x\) n*(®) 


where 


x&X 


n [»'(’ 


s(x) = 




n[»‘< 

: ,eA f 

n [n 


dx 


p{n) = 


ieAT 


n [»•(= 


ieAT 


dx 


ini*) 


n 


ieAT 


dx 


(4.14) 


(4.15) 


(4.16) 


In words, (4.14)-(4.16) amount to state that the fusion of i.i.d. cluster processes provides an 
i.i.d. cluster process whose location density s(-) is the weighted geometric mean of the agent 
location densities s*(-), while the fused cardinality p(-) is obtained by the more complicated 
expression (4.16) also involving the agent location PDFs besides the agent cardinality PMFs. 

Please notice that, in principle, both the cardinality PMF and the location PDF are 
infinite-dimensional. For implementation purposes, finite-dimensional parametrizations of 
both need to be adopted. As far as the cardinality PMF p(n) is concerned, it is enough 
to assume a sufficiently large maximum number of objects n max present in the scene and 
restrict p{ ) to the finite subset of integers {0,1,... ,n max }. As for the location PDF, two 
finitely-parameterized representations based on the SMC or, respectively, GM approaches 
are most commonly adopted. The SMC approach consists of representing location PDFs as 
linear combinations of delta Dirac functions, i.e. 


N p 

s( x ) = 5(x- Xj) (4.17) 

i=i 

Conversely, the GM approach expresses location PDFs as linear combinations of Gaussian 


73 



CHAPTER 4. DISTRIBUTED MULTI-OBJECT FILTERING 


components, i.e. 

N G 

s(x ) = J\f(x] x',, Pj) (4.18) 

3 =1 

Recent work [UJCR10] has presented an implementation of the distributed multi-object fusion 
following the SMC approach. For DMOF over a sensor network, typically characterized by 
limited processing power and energy resources of the individual nodes, it is of paramount 
importance to reduce as much as possible local (in-node) computations and inter-node data 
communication. In this respect, the GM approach promises to be more parsimonious (usually 
the number of Gaussian components involved is orders of magnitude lower than the number 
of particles required for a reasonable tracking performance) and hence preferable. For this 
reason, a GM implementation of the fusion (4.15) has been adopted in the present work 
which, further, exploits consensus in order to carry out the fusion in a fully distributed way. 
For the sake of simplicity, let us consider only two agents, a and b, with GM location densities 


n g 

s l {x) = y~)c»4 J\f(x\ x \, pf) , i = a,b 

3 =1 

A first natural question is whether the fused location PDF 



s\x) 

3 

/ [s a (*r 

s\x) 

1 — UJ 

dx 


(4.19) 


(4.20) 


is also a GM. Notice that (4.20) involves exponentiation and multiplication of GMs. To 
this end, it is useful to draw the following observations concerning elementary operations on 
Gaussian components and mixtures. 


• The power of a Gaussian component is a Gaussian component, more precisely 


[aA(i; x , P)] u = a u (3(u, P) A f(^x] x, —^ 

where 

[det(2vrP)]2 


(4.21) 

(4.22) 


• The product of Gaussian components is a Gaussian component, more precisely [WM06, 
Wil03] 

aiAf(x\xi,Pi) ■ a 2 Af(x; x 2 , T 2 ) = a 12 N(x\ x 12 , P12) (4.23) 

where 


P\ 2 

= (Pf 1 + J P 2 - 1 )' 1 

(4.24) 

X 12 

= P 1 2 (P\ X Xl + P'2 1 X2^j 

(4.25) 
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a 12 = ai a -2 - x 2 ; 0 , Pi + P 2 ) 


(4.26) 


• Due to (4.23) and the distributive property, the product of GMs is a GM. In particular, 
if s a (-) and s b (-) have Nq and, respectively, Nq Gaussian components, then s a (-) s b (-) 
will have NqNq components. 

• Exponentiation of a GM does not provide, in general, a GM. 

As a consequence of the latter observation, the fusion (4.19)-(4.20) does not provide a GM. 
Hence, in order to preserve the GM form of the location PDF throughout the computations a 
suitable approximation of the GM exponentiation has to be devised. In [J 11 IO 6 ] it is suggested 
to use the following approximation 

'n g 

^2ajjV(x-,Xj,Pj) 

J= l 

As a matter of fact, the above approximation seems reasonable whenever the cross-products 
of the different terms in the GM are negligible for all x; this, in turn, holds provided that 
the centers Xj and Xj, i / j, of the Gaussian components are well separated, as measured 
by the respective covariances Pi and Pj. In geometrical terms, the more separated are the 
confidence ellipsoids of the Gaussian components the smaller should be the error involved in 
the approximation (4.27). In mathematical terms, the conditions for the validity of (4.27) 
can be expressed in terms of Mahalanobis distance [Mah36] inequalities of the form 


n g 

- E [ a jN{x\ 

3 = 1 


^Pi)T = 


N G 

E< 

3 = 1 


/3(cu, Pj) AM x; Xj , — 


P 


(4.27) 


(Xj - Xj) T P i 1 (Xi - Xj) > 1 
(xj - Xj) T H” 1 (^ - Xj) > 1 


Provided that the use of (4.27) is preceded by a suitable merging step that fuses Gaussian com¬ 
ponents with Mahalanobis (or other type of) distance below a given threshold, the approxima¬ 
tion seems reasonable. To this end, the merging algorithm proposed by Salmond [Sal88,Sal90] 
represents a good tool. 

Exploiting (4.27), the fusion (4.20) can be approximated as follows: 


n g n g 


n g n g 


T.F < V(x; i$, i?/) ££ ag>V(x; xg>, P, 


s(x) = 


i=lj=l 

n g n g 


i=lj=l 


£2>«'W(x;xg\P j 


i=ljr = l 


dx 


n g n g 

EE- 

i=ij=i 


ab 

ij 


where 


pab 

ij 


“(f?r‘+( i-w) (f‘) 


(4.28) 


(4.29) 
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(4.30) 

(4.31) 


Notice that (4.28)-(4.31) amounts to performing a Cl fusion on any possible pair formed by a 
Gaussian component of agent a and a Gaussian component of agent b. Further the coefficient 
a“j > of the resulting (fused) component includes a factor Af(xf — x0 , P^/u + Pj / (1 — w)) 
that measures the separation of the two fusing components (x“,Pf) and ( x^Pj^j. In (4.28) 
it is reasonable to remove Gaussian components with negligible coefficients . This can be 
done either by fixing a threshold for such coefficients or by checking whether the Mahalanobis 
distance 

\j (®“ - x iV + “) ( X i - X i ) ( 4 - 32 ) 

falls below a given threshold. The fusion (4.20) can be easily extended to IV > 2 agents 
by sequentially applying the pairwise fusion (4.29)-(4.30) N — 1 times. Note that, by the 
associative and commutative properties of multiplication, the ordering of pairwise fusions is 
irrelevant. 


4.2 Consensus GM-CPHD filter 


This section presents the proposed Consensus Gaussian Mixture-CPHD (CGM-CPHD) filter 
algorithm. The sequence of operations carried out at each sampling interval k in each node 
i £ N of the network is reported in Table 4.1. All nodes i £ J\f operate in parallel at each 
sampling interval k in the same way, each starting from its own previous estimates of the 
cardinality PMF and location PDF in GM form 




( N h) k - 1 


' k ~ 1 h =i 


(4.33) 


and producing, at the end of the various steps listed in Table 4.1, its new estimates of the 
CPHD as well as of the object set, i.e. 



{(“54>q)j(!('‘ 


n 


(4.34) 


A brief description of the sequence of steps of the CGM-CPHD algorithm is in order. 


1. First, each node i performs a local GM-CPHD filter update exploiting the multi-object 
dynamics and the local measurement set Yjf. The details of the GM-CPHD update (pre¬ 
diction and correction) can be found in [VVC07]. A merging step, to be described later, 
is introduced after the local update and before the consensus phase in order to reduce 
the number of Gaussian components and, hence, alleviate both the communication and 
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the computation burden. 

2. Then, consensus takes place in each node i involving the subnetwork M l . Each node 
exchanges information (i.e., cardinality PMF and GM representation of the location 
PDF) with the neighbors; more precisely node i transmits its data to nodes j such 
that i e AG and waits until it receives data from j £ Af l \{i}. Next, node i carries 
out the GM-GCI fusion in (4.16) and (4.28)-(4.31) over A/"*. Finally, a merging step is 
applied to reduce the joint communication-computation burden for the next consensus 
step. This procedure is repeatedly applied for an appropriately chosen number L > 1 
of consensus steps. 

3. After the consensus, the resulting GM is further simplified by means of a pruning step to 
be described later. Finally, an estimate of the object set is obtained from the cardinality 
PMF and the pruned location GM via an estimate extraction step to be described later. 


Table 4.1: Distributed GM-CPHD (DGM-CPHD) filter pseudo-code 

procedure DGM-CPHD(node i, time k) 

Local Filtering 
Local GM-CPHD Prediction 
Local GM-CPHD Correction 
Merging 

Information Fusion 
for l = 1,... ,L do 

Information Exchange 
GM-GCI Fusion 
Merging 

end for 

Extraction 

Pruning 

Estimate Extraction 

end procedure 


> See Table 4.3 

> See Table 4.4 


> See (4.16) and (4.28) 
> See Table 4.2 


> See (2.70) and [VVC07] 

> See (2.71) and [VVC07] 

> See Table 4.2 


While the local GM-CPHD update (prediction and correction) is thoroughly described in 
the literature [VVC07] and the GM-GCI fusion has been thoroughly dealt with in the previous 
section, the remaining steps (merging, pruning and estimate extraction) are outlined in the 
sequel. 

Merging 

Recall [VM06,VVC07] that the prediction step (2.70) and the correction step (2.71) of the 
PHD/CPHD filter make the number of Gaussian components increase. Hence, in order 
to avoid an unbounded growth with time of such components and make the GM-(C)PHD 
filter practically implementable, suitable component reduction strategies have to be adopted. 
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In [VM06, Section III.C, Table II], a reduction procedure based on truncation of components 
with low weights, merging of similar components and pruning, has been presented. For use 
in the proposed CGM-CPHD algorithm, it is convenient to deal with merging and pruning 
in a separate way. A pseudo-code of the merging algorithm is in Table 4.2. 


Table 4.2: Merging pseudo-code 


procedure Merging({x;, P t , a j}^ , j m ) 
t = 0 

1 = {1,...,N G } 
repeat 
t = t+ 1 
j = argrnax a* 


MC = { 


el 


(Xj 



<*t= J2 


OLi 


i£MC 

x t = — V ctiXi 

r\: i z ' 


1 


ieMC 


k ieMC 

1 = 1\MC 

until 1 0 

return 


a,- 


Pi + (Xk - Xi) (x k - Xi) 


end procedure 




Pruning 

To prevent the number of Gaussian components from exceeding a maximum allowable value, 
say N max , only the first N max Gaussian components are kept while the other are removed 
from the GM. The pseudo-code of the pruning algorithm is in Table 4.3. 

Table 4.3: Pruning pseudo-code 


procedure PRUNiNG({xj, Pi, a*}^ , N max ) 
if N g > N max then 

1 = {indices of the N max Gaussian components with highest weights cti} 
return {x*, Pi, Oj} i6X 

end if 

return {x,;, P t , ai}^l{ ax 

end procedure 


Since the CPHD filter directly estimates the cardinality distribution, the estimated num- 
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ber of objects can be obtained via Maximum A Posteriori (MAP) estimation, i.e. 

fik = maxpfc(n) (4.35) 

Given the MAP-estimated number of objects hk, estimate extraction is performed via the 
algorithm in Table 4.4 [VM06, Section III.C, Table III]. The algorithm extracts the hk local 
maxima (peaks) of the estimated location PDF keeping those for which the corresponding 
weights are above a preset threshold r y e . 


Table 4.4: Estimate extraction pseudo-code 

procedure Estimate Extraction^®*, Pi, ,n/c,7 e ) 

X = {indices of the hk Gaussian components with highest weights 

x t = <D 

for i G X do 

if aihk > 7e then 
X t = X t U Xi 

end if 
end for 
return Xt 
end procedure 


4.3 Performance evaluation 


To assess performance of the proposed CGM-CPHD algorithm described in section 4.2, a 
2-dimensional (planar) multi-object tracking scenario is considered over a surveillance area 
of 50 x 50 [km 2 ], wherein the sensor network of Fig. 4.1 is deployed. The scenario consists of 
6 objects as depicted in Fig. 4.2. 

The object state is denoted by x = \p x , p x , p y , p y ] T where ( p x ,p y ) and (PxLPy) represent 
the object Cartesian position and, respectively, velocity components. The motion of objects 
is modeled by the filters according to the nearly-constant velocity model: 


1 
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rn 2 
■ L s 


(4.36) 


where a w = 2 [m/s 2 ] and the sampling interval is T s = 5 [s]. 

As it can be seen from Fig. 4.1, the sensor network considered in the simulation consists 
of 4range-only ( Time Of Arrival, TOA) and 3 bearing-only ( Direction Of Arrival, DOA) 
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xIO 4 



Figure 4.1: Network with 7 sensors: 4 TOA and 3 DOA. 


sensors characterized by the following measurement functions: 

[ Z[(p x -x l ) +j (p y -y% DOA 
ti(x) = l __ (4.37) 

[ v {px - X 1 ) 2 + (p y - y*) 2 , TOA 

where (x l ,y z ) represents the known position of sensor i. The standard deviation of DOA 
and TOA measurement noises are taken respectively as ctdoa = 1 [°] and ®toa = 100 [m\. 
Because of the non linearity of the aforementioned sensors, the UKF [JU04] is exploited in 
each sensor in order to update means and covariances of the Gaussian components. 

Clutter is modeled as a Poisson process with parameter A c = 5 and uniform spatial 
distribution over the surveillance area; the probability of object detection is Pd = 0.99. 

In the considered scenario, objects pass through the surveillance area with no prior infor¬ 
mation for object birth locations. Accordingly, a 40-component GM 

40 

db(x) = a j N{x] Xj , Pj) (4.38) 

3 = 1 

has been hypothesized for the birth intensity. Fig. 4.3 gives a pictorial view of c4(x); notice 
that the center of the Gaussian components are regularly placed along the border of the 


80 










4.3. PERFORMANCE EVALUATION 


5 


4 


3 


2 


1 


0 


Figure 4.2: object trajectories considered in the simulation experiment. The start/end point 
for each trajectory is denoted, respectively, by «\B. 

surveillance region. For all components, the same covariance Pj = diag(10 6 , 10 4 ,10 6 , 10 4 ) 
and the same coefficient olj = a, such that rib = 40a is the expected number of new-born 
objects, have been assumed. 

The proposed consensus CGM-CPHD filter is compared to an analogous filter, called 
Global GM-CPHD (GGM-CPHD), that at each sampling interval performs a global fusion 
among all network nodes. Multi-object tracking performance is evaluated in terms of the Op¬ 
timal SubPattern Assignment (OSPA) metric [SVV08]. The reported metric is averaged over 
N mc = 200 Monte Carlo trials for the same object trajectories but different, independently 
generated, clutter and measurement noise realizations. The duration of each simulation trial 
is fixed to 500 [s] (100 samples). 

The parameters of the CGM-CPHD and GGM-CPHD filters have been chosen as follows: 
the survival probability is Ps = 0.99; the maximum number of Gaussian components is 
N max = 25; the merging threshold is 7 m = 4; the truncation threshold is 7 1 = 10 , the 

extraction threshold is 7 e = 0.5; the weight of each Gaussian component of the birth PHD 
function is chosen as aj = 1.5 • 10~ 3 . 

Figs. 4.4, 4.5, 4.6 and 4.7 display the statistics (mean and standard deviation) of the 
estimated number of objects obtained with CGM-CPHD (Fig. 4.4) and CGM-CPHD with 
L = 1 (Fig. 4.5), L = 2 (Fig. 4.6) and L = 3 (Fig. 4.7) consensus steps. Fig. 4.8 reports 
the OSPA metric (with Euclidean distance, p = 2, and cutoff parameter c = 600) for the 
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M 


xIO 4 


Figure 4.3: Borderline position GM initialization. The symbol o denotes the component mean 
while the blue solid line and the red dashed line are, respectively, their 3a and 2 a confidence 
regions. 


same filters. As it can be seen from Fig. 4.8, the performance obtained with three consensus 
steps is significantly better than with a single one, and comparable with the one given by 
the non scalable GGM-CPHD filter which performs a global fusion over all network nodes. 
Similar considerations hold for cardinality estimation. These results show that by applying 
consensus, for a suitable number of steps, performance of distributed scalable algorithms is 
comparable to the one provided by the non scalable GGM-CPHD. 

To provide an insightful view of the consensus process, Fig. 4.9 shows the GM repre¬ 
sentation within a particular sensor node (TOA sensor 1 in Fig. 4.1) before consensus and 
after L = 1, 2, 3 consensus steps. It can be noticed that consensus steps provide a progressive 
refinement of multi-object information. In fact, before consensus the CPHD reveals the pres¬ 
ence of several false objects (see Fig. 4.9a); this is clearly due to the fact that the considered 
TOA sensor does not guarantee observability. Then, in the subsequent consensus steps (see 
Figs. 4.9b-4.9d), the weights of the Gaussian components corresponding to false objects as 
well as the covariance of the component relative to the true object are progressively reduced. 
Summing up, just one consensus step leads to acceptable performance but, in this case, the 
distributed algorithm tends to be less responsive and accurate in estimating object number 
and locations when compared to the centralized one (GGM-CPHD). However, by performing 
additional consensus steps, it is possible to retrieve performance comparable to the one of 
GGM-CPHD. 
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Figure 4.4: Cardinality statistics for GGM-CPHD. 
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Figure 4.5: Cardinality statistics for CGM-CPHD with L = 1 consensus step. 


In order to have an idea about the data communication requirements of the proposed 
algorithm, it is useful to examine the statistics of the number of Gaussian components of 
each node just before carrying out a consensus step; this determines, in fact, the number of 
data to be exchanged between nodes. To this end, Table 4.5 reports the average, standard 
deviation, minimum and maximum number of components in each node and at each consensus 
step during the simulation. Notice that, in the considered case-study, the average number of 
Gaussian components to be transmitted is about 20 -f- 30. Since each component involves 15 
real numbers (4 for the mean, 10 for the covariance and 1 for the weight), the overall average 
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Figure 4.6: Cardinality statistics for CGM-CPHD with L = 2 consensus steps. 



communication load per consensus step is about 1.2 -4- 1.8 kbytes if a 4-byte single-precision 
floating-point representation is adopted for each real number. 

To give a rough idea of the computation time, our MATLAB implementation on a Personal 
Computer with 4.2 GHz clock exhibited an average processing time per consensus step in 
the range of 45 4- 55 ms. Clearly, this time can be significantly reduced with a different, e.g. 
C/C++, implementation. 

The CCPHD filter has also been successfully applied in scenarios with range and/or 
Doppler sensors [BCF + 13b] and with passive multi-receiver radar systems [BCF + 13c]. 
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Figure 4.8: Performance comparison, using OSPA, between GGM-CPHD and CGM-CPHD 
respectively with L = 1, L = 2 and L = 3 consensus steps. 


Table 4.5: Number of Gaussian components before consensus for CGM-CPHD 
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AVG 
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DOAi 
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Figure 4.9: GM representation within TOA sensor 1 (see Fig. 2) at the initial time instant 
throughout the consensus process. Means and 99.7% confidence ellipses of the Gaussian 
components of the location PDF are displayed in the horizontal plane while weights of such 
components are on the vertical axis. 
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In this chapter a new approximation of the d-GLMB filter is presented [FVPV15]. It 
will be shown that the GLMB distribution can be used to construct a principled approx¬ 
imation to an arbitrary labeled RFS density that matches the PHD and the cardinality 
distribution [PVV + 14], The resulting filter is referred to as the Marginalized d-GLMB (M<5- 
GLMB) tracking filter since it can be interpreted as a marginalization over the data asso¬ 
ciations. The proposed filter is consequently computationally cheaper than the (5-GLMB 
filter while still preserving key summary statistics of the multi-object posterior. Impor¬ 
tantly, the M5-GLMB filter facilitates tractable multi-sensor multi-object tracking. Unlike 
PHD [VM06], CPHD [VVC07] and Multi-Bernoulli based filters [VVC09], the proposed ap¬ 
proximation accommodates statistical dependence between objects. An alternative derivation 
of the LMB filter proposed in [RVVD14] based on the newly proposed Md-GLMB filter is 
presented [FVPV15]. 
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5.1 The (5-GLMB filter 

An efficient approach to multi-object tracking was presented in [VV13] using the 5-GLMB 
distribution (2.86) (or equivalently (2.87)), i.e. 


TT(X) = A(X) 

E ^ (/ ’ ?) 5/(£(X)) [p«)] X , 

(5.1) 


(7,0e.F( L)xS 


= A(X) 


(5.2) 


/eJ’(L) £es 


The 5-GLMB density naturally arises in multi-object tracking problems when using the 
standard detection based measurement model . In the following, the prediction and update 
steps of the 5-GLMB filter is briefly recalled; details can be found in [VV13, VVP14]. 

5.1.1 5-GLMB prediction 

The standard multi-object dynamic model is described as follows. Given the multi-object 
state Z, each state (£,£-) € Z either continues to exist at the next time step with proba¬ 
bility Pg(£,£-) and evolves to a new state ( x ,£) with probability density Tk\k-i( x \Ci ■£-), or 
disappears with probability 1 — Ps((,£-)- Note that the label of the objects is preserved in 
the transition, only the kinematic part of state changes. Assuming that Z has distinct labels 
and that conditional on Z, the transition of the kinematic states are mutually independent, 
then the set W of surviving objects at the next time is a labeled multi-Bernoulli RFS [VV13] 

/s(W|Z) = A(W) A(Z) l £(z) (£(W)) [<1 (W; -)] Z , (5.3) 

where 

T(W;C,G)= ^ Se-W Ps ((,£-) f (x\(,£-)+[l-1 C(W) (£-)\(1-P S ((,£-)) . (5.4) 

The A(Z) in (5.3) ensures that only Z with distinct labels are considered. 

The set of new objects born at the next time step is distributed according to 

/b(Y) = A(Y) wb{E(Y)) [pb] Y (5.5) 

The birth density f B {-) is defined onXxL and f B ( Y) = 0 if Y contains any element x with 
£(x) ^ L. The birth model (5.5) includes both labeled Poisson and labeled multi-Bernoulli 
densities. The multi-object state at the next time X is the superposition of surviving objects 
and new born objects, i.e. X = W U Y. Since the label spaces L and B are disjoint, the 
labeled birth objects and surviving objects are independent. Thus the multi-object transition 
density turns out to be the product of the transition density (5.3) and the density of new 
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objects (5.5) 

/(X|Z) = fs(X n (X x L) |Z) f B (X - (X x L)) . (5.6) 

Additional details can be found in [VV13, Subsection IV.D]. 

If the current multi-object prior density is a 5-GLMB of the form (2.86), then the multi¬ 
object prediction density is also a d-GLMB given by 


7T 


fc|fc-i(X) = A(X) E ^//(AX)) 

(/,C)eJ(L)x= 


Pk\l-i 


where 


(5.7) 


w 


V*) _ 

k\k 


^ flL_) , 


Pk\k- iOM) = _(£)ps\x,£) + 1 m{£)pb{x,£) , 

P S OM) = ~ - 


Vs\t) 

v < £ ) W = (psM,p$1 1 m) , 

4P(L) = [ V f] L E lj(L)[l-v ( s°} J - L wi J ^ 


J CL_ 


(5.8) 

(5.9) 

(5.10) 

(5.11) 

(5.12) 


5.1.2 5-GLMB update 

The standard multi-object observation model is described as follows. For a given multi-object 
state X, each state x 6 X is either detected with probability Po (x) and generates a point 
y with likelihood g(y |x), or missed with probability 1 — Po (x), i.e. x generates a Bernoulli 
RFS with parameter (P^(x), g(-|x)). Assuming that conditional on X these Bernoulli RFSs 
are independent, then the set IF C Y of detected points (non-clutter measurements) is a 
multi-Bernoulli RFS with parameter set {(Pd(x), <?(• |x)) : x £ X}. The set C C Y of false 
observations (or clutter), assumed independent of the detected points, is modeled by a Poisson 
RFS with intensity function «(•). The multi-object observation Y is the superposition of the 
detected points and false observations, i.e. Y = W U C. Assuming that, conditional on X, 
detections are independent, and that clutter is independent of the detections, the multi-object 
likelihood is given by 


g k {Y |X) = E [Vy(-;0)] X , (5.13) 

6»e0(£(X)) 

where 0(/) is the set of mappings 6 : I —> {0,1, ...,M}, such that 9{i) = 6{i') > 0 implies 
i = i' , and 

r p D (x,e) 9 t (vm\x,e) | ite{t)>0 

ip Y {x,P,e) = < K {ye(i)) . (5.14) 

(l -P D {x,i), if Q(t) = 0 
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Note that an association map 6 specifies which tracks generated which measurements, i.e. 
track I generates measurement ygm G Y, with undetected tracks assigned to 0. The condition 
“9(i) = 9{i') > 0 implies i = i h \ means that a track can generate at most one measurement, 
and a measurement can be assigned to at most one track, at one time instant. Additional 
details can be found in [VV13, Subsection IV.C]. 

If the current multi-object prediction density is also a 5-GLMB of the form (2.86), then 
the multi-object posterior density is a d-GLMB given by 


w k (X) = A(X) 


E E 

(GOe-FOLjxS 0e©(/) 






lX 


(5.15) 


where 0(1) denotes the subset of the current maps with domain I, and 


w k^ d \ Y k) « 

= 

Pk ,8) (M) = 


5 

(5.16) 

(5.17) 

(5.18) 

Vyf\ £ ) 



( P D {x,£) g{y m \x,£) 

I / \ 5 

if 9(£) > 0 



l K(jfe(/)) 

(5.19) 


(l -P D (x,i), 

O 

II 

Set 



Notice that the new association maps 0 can be added (stacked) to their respective association 
histories £ in order to have again the more compact form (2.86) for the updated <5-GLMB 
(5.15). 


5.2 GLMB approximation of multi-object densities 

In this section a GLMB approximation of the multi-object density with statistically dependent 
objects is proposed. In particular a GLMB density that matches the multi-object density of 
interest in both the PHD and cardinality distribution is derived. The strategy is inspired by 
Mahler’s iid cluster approximation in the CPHD filter [Mah07a], which has proven to be very 
effective in practical multi-object filtering problems [Mahl4,GSW09,SWS09]. Moreover, the 
GLMB approximation captures object dependencies whereas the iid cluster approximation 
assumes that individual objects are iid from a common single-object density. Proof of the 
result is given in the appendix A. 

Our result follows from the observation that any labeled RFS density 7r(-) on ^(XxL) 
can be written as 

7r(X) = to(£(X))p(X) (5.20) 


90 



5.3. MARGINALIZATIONS OF THE 6-GLMB DENSITY 


where 


w({£i,... ,4}) = j ..., (x n ,£ n )})dx i • • • 

P{{(X l,£l),---,(x n ,£n)}) = 


dx r 


a 4), ■••,(*«> 4)}) 


u>({£ i,. • .,4}) 


(5.21) 

(5.22) 


It is implicitly assumed that p(X) is defined to be zero whenever w(£(X.)) is zero. 

Note that since 7r(-) is symmetric in its arguments, the integral in (5.21) and w is symmet¬ 
ric in £\,... ,4- Moreover, since Y^LeT( L) w(L) = 1, w is indeed a probability distribution on 
J^L) and can be interpreted as the probability that the labeled RFS has label set { £\ ,..., 4 }- 
For w{{£\,... , 4 }) > 0, p({(xi,4), • • •, (x' n ,4)}) is the joint probability density (on X n ) of 
the kinematic states xi, ... ,x n given that their corresponding labels are £i,... ,£ n . 


Proposition 1 (GLMB approximation). 

The GLMB that matches the cardinality distribution and the PHD of a labeled RFS with 
density 7r(-) is given by: 


7f(X) = 


E ™ (/} (AX)) 


i( 7 ) 


(5.23) 


I&F{ L) 


where, 


w^\L) = 5 I (L)w(I) (5.24) 

p (I \x,£) = l/(£)p AW (x,f) (5.25) 

V{h,...,t n }{x,£) = J p({(x,£), (xi,£i),..., (x n ,£ n )})d(x l ,.. .,x n ) (5.26) 

Remark 13. Note that in [BVV14, Section V] a GLMB was proposed to approximate a 
particular family of labeled RFS densities that arises from multi-object filtering with merged 
measurements. By applying Proposition 1, it can be shown that the approximation used 
in [BVV14, Section V] preserves the cardinality distribution and first moment. 

Remark 14. In multi-object tracking, the matching of the labeled PHDs d{-,£) and d(-,£) 
in Proposition 1 is a stronger result than simply matching the unlabeled PHDs alone. 


5.3 Marginalizations of the 5-GLMB density 

One of the main factors increasing the computational burden of the d-GLMB filter [VVP14] 
is the exponential growth of the number of hypotheses in the update of the prior (5.15), which 
gives rise to the explicit sum over an association history variable £. Further, note that the 
number of association histories is considerably increased in a multi-sensor context. The idea 
behind the Md-GLMB filter is to construct a principled GLMB approximation 7r(-) to the 
posterior density 7r(-) by means of a marginalization with respect to the association histories 
£. In this section, the time index is dropped for the sake of simplicity and the attention is 
devoted to the Mci-GLMB approximation corresponding to any (5-GLMB density. 
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Proposition 2 (The M5-GLMB density). 

A Marginalized 5-GLMB (M5-GLMB) density J k{ ) corresponding to the 5-GLMB density 
7r(*) in (2.86) is a probability density of the form 


tt(X) = A(X) £ 5/(£(X)) W « 

ieT( L ) 

£GS 



p (I \x,£) = —tty . 


(5.27) 

(5.28) 

(5.29) 


TTie M5-GLMB (5.27) preserves both (labeled) PHD and cardinality distribution of the orig¬ 
inal 5-GLMB density 7r(-) in (2.86). 


Note now that the M5-GLMB density provided by Proposition 2 is a 5-GLMB which 
aggregates the association histories £ relative to the label set I via (5.28)-(5.29). In fact, the 
marginalized density (5.27) has no association history set H. 

The M5-GLMB density can be exploited to construct an efficient recursive multi-object 
tracking filter by calculating the M5-GLMB approximation step after the 5-GLMB update, 
and predicting forward in time using the 5-GLMB prediction. 


5.4 The M5-GLMB filter 

The M5-GLMB filter propagates a M5-GLMB multi-object posterior density forward in time 
via the multi-object Bayesian recursion (2.97)-(2.98). Since the M5-GLMB is still a 5-GLMB, 
the equations derived for the 5-GLMB filter [VVP14] as a closed-form solution for the eqs. 
(2.97)-(2.98) still hold. In the following, the prediction and update steps for the 5-GLMB 
filter are briefly recalled; additional details can be found in [VVP14,FVPV15]. 

5.4.1 M5-GLMB prediction 

Given the multi-object state Z, each state ((,£-) € Z either continues to exist at the next 
time step with probability Ps{(, £-) and evolves to a new state (x, £) with probability density 
Pklk- i(ic|C, ■£-)> or dies with probability 1 — Ps((, £-)■ The set of newborn objects at the next 
time step is distributed according to the LMB 

fB={( S S,p'S)} le , (5-30) 

with L_nB = 0 and tcs(X) defined in (2.91). The multi-object state at the current time X^ 
is the superposition of surviving objects and new born objects. 

If the current multi-object prior density is a M5-GLMB of the form (5.27), then the 
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multi-object prediction density is also a M5-GLMB given by 


7T fc | fe —,(X) = A(X) E S T (£(X)) 

ieT( L) 


w 


(i) 

k\k-l 


Pk\l-1 



where 


w k\l-l = ^b(A l -) 4 \ I n L -) > 

Pk\l-i( x A) = !l_ + 1b (£)pb{x,£) , 

4 AM) =-(Tw-’ 

A g) 

4 7) w = (a(-a),p2 1 (a)) , 

4 7) (A = [4 J) ] L E a(A[i-4 7) ] j - l 4 j) 5 

JCL_ 


(5.31) 


(5.32) 

(5.33) 

(5.34) 

(5.35) 

(5.36) 


which is exactly the 5-GLMB prediction step (5.7) with no association histories from previous 
time step, i.e. H = 0, and with the convention of having the superscript (!) instead of (£) 
due to the marginalization (5.28)-(5.29). 

Remark 15. The number of components (4+ 1 |fc , Pfe+i|fc) corn P u ted after the Md-GLMB 
prediction step (5.31) is |J r (Lo : fc+i)|- On the other hand, the number of components (4+i]fc; 
Pfe+i| fe ) after the d-GLMB prediction (5.7) is ^(Lo^+i) x S| for and l s l for Pk+i\k- 

Notice that the number of weights 4+i|fc the M(5-GLMB is significantly lower than the 
4+p k 5-GLMB. As for the number of location PDFs p\J +l | fc , it is worth noticing that the 
growth rate of the association histories £ £ H is super-exponential with time [VV13, VVP14], 
while the growth rate of the cardinality of J-^Lo^+i) is exponential. 

The use of the M5-GLMB approximation further reduces the number of hypotheses in 
the posterior density while preserving the PHD and cardinality distribution. Moreover, the 
M5-GLMB is in a form that is suitable for efficient and tractable information fusion (i.e. 
multi-sensor processing) as will be shown in the next subsection. 


5.4.2 M5-GLMB update 

Given a multi-object state X, each state (. x , £) G X is either detected with probability Pd (x, £) 
and generates a measurement y with likelihood gk(y\x,£), or missed with probability 1 — 
Pd(x, £)■ The multi-object observation Y = {y\ ,..., ijm\ is the superposition of the detected 
points and Poisson clutter with intensity function «(•). 

If the current multi-object prediction density is a M5-GLMB of the form (5.27), then the 
multi-object posterior density is a 5-GLMB given by 


7T fc (X) = A(X) E 


E 


5i{C{X))w^ e \Y k ) 


Pk 


ieT( L) eee(i) 


(5.37) 
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where 0(1) denotes the subset of the current maps with domain I, and 


w k’ 0 ) (Y k ) oc wjjl-i. [vy^W ] 1 , 

4^(0 = {Pk\ k -i(-,V)AY k (-A\0)) 

<i,o) t , Pk\l-i( x ^)^Y k (x,P,e) 

p k (*) =- myy, -’ 

VY k U) 

{ Pp(x,£)g(y e ^)\x,£) 
K (ye(i)) 

1 ~Pd(x,£), 


if 6(£) > 0 
if 6(£) = 0 


(5.38) 

(5.39) 

(5.40) 

(5.41) 


Using (5.28)-(5.29), the Md-GLMB density corresponding to the <5-GLMB density in (5.37) 
is a probability density of the form (5.27) with 

U’= E «T ) > (5.42) 

0e0(J) 

Pk\ x > £ ) = i/W-Tn w k’ 0 ) Pk' e) ( x ( 5 - 43 ) 

w k eee(i) 

The M<5-GLMB density provided by (5.42)-(5.43) preserves both PHD and cardinality distri¬ 
bution of the original <5-GLMB density. 

Remark 16. Each hypothesis I £ J- (L) generates a set of |@(/)| new measurement-to-track 
association maps for the 5-GLMB posterior. The number of components (444’ 4+7|fc) 
stored/computed after the M5-GLMB update step (5.37) is |.F(]L)| • E/e-FQL) |©(-0|- On the 
other hand, the number of hypotheses (44|fc -P'k+\\k) a fter the <5-GLMB update (5.15) is 

(^(L) x H|-E/g^(L) |6(/)| for w k+i\k and l S l'E/e^(L) l©G0l ioT PkS\ k - The same conclusions 
along the lines of Remark 15 hold. 

Remark 17. After the marginalization procedure (5.42)-(5.43) only |.F(L)| hypotheses are 
retained, as all the new contributions provided by the association maps |@(/)| are aggregated 
in a single component. Notice that |.F(]L)| is the same number of hypotheses produced during 
the prediction step (5.31) (see Remark 15). Thus, the prediction step (5.37) sets the upper 
bound of the hypotheses that will be retained after each full M5-GLMB step. 

In the multi-sensor setting, for each time instant k, the update step (5.37) (or (5.15) for 
the d-GLMB) is repeatedly evaluated using different measurement sets Y l provided by the 
sensors i £ J\f (see eq. (2.101)). From Remarks 16 and 17, Md-GLMB is preferable to 5- 
GLMB in terms of memory and computational requirements, since the number of remaining 
hypotheses after each sensor update step in (2.101) is always set to |J r (IL)|. Note that this 
does not apply to d-GLMB due to the super-exponential memory/computational growth 
as reported in Remark 15. This is an important property of M5-GLMB since it yields a 
principled approximation which greatly decreases the need of pruning hypotheses with respect 
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to (5-GLMB [VVP14], In fact, pruning in the <5-GLMB might lead to poor performance in 
multi-sensor scenarios with low SNR (e.g. high clutter intensity, low detection probability, 
etc.) and limited storage/computational capabilities. For instance, this may happen if a 
subset of the sensors do not detect one or multiple objects and hypotheses associated to 
the true tracks are removed due to pruning. Furthermore, from a mathematical viewpoint, 
pruning between corrections generally produces a less informative and order-independent 
approximation of the posterior distribution in (2.101). 

5.5 The LMB filter 

The LMB filter introduced in [RVVD14] is a single component approximation of a GLMB 
density (2.77) that matches the unlabeled PHD. An alternative derivation of the LMB ap¬ 
proximation (2.90) first proposed in [RVVD14] through a connection with the MAGLMB 
approximation is provided hereafter. 

5.5.1 Alternative derivation of the LMB Filter 

Recall that a LMB density is uniquely parameterized by a set of existence probabilities 
and corresponding track densities p^\-): 

7r(X) = A(X) w(£(X))p x (5.44) 

where 

ll ,(i)=n(i-r w )n 1 rw <5-45) 

isl eeL 

p(x,£) = p( e \x) (5.46) 

By extracting individual tracks from the Md-GLMB approximation 7r(-) (5.27) it is pos¬ 
sible to prove that the same expressions for the existence probabilities and state densities 
originally proposed for the LMB filter in [RVVD14] can be obtained. In fact 


W = PiA(f G £(X)) 

(5.47) 

= ™( L ) 

(5.48) 

t&L 


= 1 dt) w ( L ) 

(5.49) 

ieT( L) 



(5.50) 

/e-F(L) 


= 1 , 

(5.51) 

(J,0£F(L)xS 
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pW(x) = 


Pr*(0M)eX) 

Pi-*(£ € £(X)) 

(5.52) 

f n((x,£)U X)5X 


J 

y(L) 

(5.53) 

d(x , £) 
r (t) ’ 

(5.54) 


where the notation for the numerator in (5.52) is defined as per [Mah07b, eq. 11.111], while 
the numerator of (5.54) follows from [Mah07b, eq. 11.112]. Notice that the numerator is 
precisely the PHD d(-) corresponding to tt(-), which by Proposition 1 exactly matches the 
PHD d(-) corresponding to 7r(-). It can be easily verified that 


d(x,£) = ^2 1 i(£)w^p^\x,£) (5.55) 

l£T{ L) 


and consequently 


N 7 y 1 i(£)w {I] p [I \x,£) 

r /e-F(L) 

(5.56) 

~(L) J2 1 i{£)w^^p^\x,£) 

r (;,5)eJ(L)x = 

(5.57) 


Remark 18. The property of matching the labeled PHD of the <5-GLMB does not hold for the 
LMB tracking filter, as shown in [RVVD14, Section III], due to the imposed multi-Bernoulli 
structure for the cardinality distribution. This can be also proved by combining (5.44) with 
Proposition 1, thus showing again that only the unlabeled PHD of the LMB matches exactly 
the unlabeled PHD of the original d-GLMB density. 


As suggested by its name, the LMB filter propagates an LMB multi-object posterior den¬ 
sity forward in time [RVVD14]. Its recursion will be briefly reviewed in the next subsections. 
Additional details on its implementation can be found in [RVVD14]. 


5.5.2 LMB prediction 

Given the previous multi-object state Z, each state (£,!'-) e Z either continues to exist at the 
next time step with probability Ps((,£~) and evolves to a new state (x,£) with probability 
density (p k \ k _ 1 (x\(,£-), or dies with probability 1 — Ps((,£-)- The set of new objects born 
at the next time step is distributed according to 

/^{(-M)}^ (5-58) 

with L_nB = 0. The multi-object state at the current time X& is the superposition of 
surviving objects and new born objects. 

Suppose that the multi-object posterior density at time A: — 1 is LMB with parameter 
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set 7Tfc_i = ^ . Then, it has been shown in [RVVD14] that the multi-object 

predicted density is also LMB with parameter set 

7r ‘l‘-i = {(4 ' ) .Ps ) )} kl _ u {(r™,p“)} KB , (5.59) 

where 

r s" 1 = Vs(£) r ^ ) (5.60) 

(7), s (^Ps('i£)Pk\k-l{ x \',fyiP('r£)') 

py(x) = - — -, (5.61) 

r, s (£) = (P s (;£),p(;£)) • (5-62) 


5.5.3 LMB update 

Given a multi-object state X, each state (x, £) € X is either detected with probability Pd (x, £) 
and generates a measurement y with likelihood gk(y\x,£), or missed with probability 1 — 
Pd(x,£). The multi-object observation Y = {j/i, ...,2 /m} is the superposition of the detected 
points and Poisson clutter with intensity function «(•). Assuming that, conditional on X, 
detections are independent, and that clutter is independent of the detections, the multi-object 
likelihood is the same of the M5-GLMB Update in (5.14). 

Suppose that the multi-object predicted density is LMB with parameter set 

*k\k-i = { ( r k\k-vPk\k-i) } teL • ( 5 - 63 ) 

The multi-object posterior density cannot be evaluated in a closed form [RVVD14], Specif¬ 
ically, the LMB 7r^,|fc_ 1 has to be converted to a 5-GLMB density by (possibly) generating 
all possible combinations of label sets I G <F(£(X)) from the constituent Bernoulli com¬ 
ponents [RVVD14]. Further, the 5-GLMB 7r fc | fc _ 1 (-) resulting from the conversion process 
undergoes the <j-GLMB update step (5.15). Finally, the GLMB posterior 7Tfc(-) is approxi¬ 
mated by an LMB that matches exactly the first moment of the unlabeled posterior intensity. 
Thus, the updated LMB is given by 


where 


7r l (X) = {(r<' ) ,pf)} t£L . (5.64) 

4 <> = £ 1,(<)«L" ) (U), (5.65) 

(/,9)eJ(L)x0(;) 

Pk\ x ) = -Je) Y 1 i( £ ) w k ,e \ Y k)Pk ) (M), (5-66) 

r k (/,0)eJ’(L)x0(7) 

w [I ’ e) (Y k ) oc [rff}] 1 1] i 1 - r k\k-i) II , (5-67) 

7eL\7 iei 
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IS), Pk\k-i(x,t)i>Y k (x,P,9) 

Pk *) =-- 




ipY k (x,e;8) 


! Pp(x,£)g(y ew \x,£) 
«(i !e{i)) 

1 ~Pd(x,£), 


if 9(£) > 0 
if 6{£) = 0 


(5.68) 

(5.69) 

(5.70) 


Additional details for an efficient implementation of the LMB filter can be found in [RVVD14]. 


5.6 Performance evaluation 

To assess performance of the proposed Marginalized 5-GLMB (Md-GLMB) tracker, a 2- 
dimensional multi-object tracking scenario is considered over a surveillance area of 50 x 
50 [km 2 ] . Two sensor sets are used to represent scenarios with different observability capa¬ 
bilities. In particular: 

i. a single radar in the middle of the surveillance region is used as it guarantee observ¬ 
ability; 

ii. a set of 3 range-only (Time Of Arrival, TOA) sensors, deployed as shown in Fig. 5.1, are 
used as they do not guarantee observability individually, but information from different 
sensors need to be combined to achieve it. 


The scenario consists of 5 objects as depicted in Fig. 5.2. For the sake of comparison, 
the M5-GLMB is also compared with the d-GLMB (d-GLMB) [VV13, VVP14] and LMB 
(LMB) [RVVD14] filters. The three tracking filters are implemented using Gaussian Mixtures 
to represent their predicted and updated densities [VVP14,RVVD14]. Due to the non linearity 
of the sensors, the Unscented Kalman Filter (UKF) [JU04] is exploited to update means and 
covariances of the Gaussian components. 

The kinematic object state is denoted by x = \p x , p x , p y , p y \ T , he. the planar position 
and velocity. The object motion is modeled according to the Nearly-Constant Velocity (NCV) 
model [FS85,FS86, BSF88, BSL95]: 


1 
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CO 


(5.71) 


where: Wk is a white noise with zero mean and covariance matrix Q, c t w = 5 [m/s 2 ] and the 
sampling interval is T s = 5 [s]. 


98 



5.6. PERFORMANCE EVALUATION 



— Surv. Area # TOA sensor 


Figure 5.1: Network with 3 TOA sensors. 


The radar has the following measurement function: 


h(x) = 


Z[(p x - x r ) + j (p y - y r )] 
V (Px - x r ) 2 + ( Py - y r ) 2 


(5.72) 


where ( x r ,y r ) represents the known position of the radar and the standard deviation of 
the measurement noise is erg = 1 [°] for the azimuth and oy = 100 [rn] for the range. The 
measurement functions of the 3 TOA of Fig. 5.1 are: 


h(x) = \J (; p x - X s ) 2 + (p y - y s Y 


(5.73) 


where ( x s ,y s ) represents the known position of sensor (indexed with) s. The standard devi¬ 
ation of the TOA measurement noise is taken as cttoa = 100 \m]. 

The clutter is characterized by a Poisson process with parameter A c = 15. The probability 
of object detection is Pd = 0.85. 

In the considered scenario, objects pass through the surveillance area with partial prior 
information for object birth locations. Accordingly, a 10-conrponent LMB RFS ttb = 

^ has been hypothesized for the birth process. Table 5.1 gives a detailed summary 
of such components. Due to the partial prior information on the object birth locations, 
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xIO 4 



Figure 5.2: Object trajectories considered in the simulation experiment. The start/end point 
for each trajectory is denoted, respectively, by «\B. The -k indicates a rendezvous point. 


Table 5.1: Components of the LMB RFS birth process at a given time k. 

r-M = 0.09 

Pb\ x ) = m^\P B ) 

P B = diag(l0 6 ,10 4 ,10 6 ,10 4 ) 


Label 

(k, 1) 

(k, 2) 

(k, 3) 

(£) 

m B 

[0, 0, 40000, 0] T 

[0, 0, 25000, 0] T 

[0, 0, 5000, 0] T 

Label 

(k, 4) 

(k, 5) 

(k, 6) 

(£) 

m B 

[5000, 0, 0, 0] T 

[25000, 0, 0, 0] T 

[36000, 0, 0, 0] T 


Label 

(*, 7) 

(k, 8) 

{t) 

m B 

[50000, 0, 15000, 0] T 

[50000, 0, 40000, 0] T 


Label 

(k, 9) 

(k, 10) 

(£) 

m B 

[40000, 0, 50000, 0] T 

[10000, 0, 50000, 0] T 


some of the LMB components cover a state space region where there is no birth. Clutter 
measurements are, therefore, more prone to generate false objects. 
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Multi-object tracking performance is evaluated in terms of the Optimal SubPattern Anal¬ 
ysis (OSPA) metric [SVV08] with Euclidean distance p = 2 and cutoff c = 600 [rn]. The 
reported metric is averaged over 100 Monte Carlo trials for the same object trajectories but 
different, independently generated, clutter and measurement noise realizations. The duration 
of each simulation trial is fixed to 1000 [s] (200 samples). 

The three tracking filters are coupled with the parallel CP HD look ahead strategy described 
in [VV13, VVP14], The CPHD [VVC07] filter. 

5.6.1 Scenario 1: radar 

Figs. 5.3, 5.4 and 5.5 display the statistics (mean and standard deviation) of the estimated 
number of objects obtained, respectively, with M5-GLMB, 5-GLMB and LMB tracking filters. 
As it can be seen, all the algorithms estimate the object cardinality accurately, with no 
substantial differences. This result implies that, in the presence of a single sensor guaranteeing 
observability, the approximations made by both M5-GLMB and LMB trackers are not critical 
in that they provide performance comparable to 5-GLMB with the advantage of a cheaper 
computational burden and reduced storage requirements. Note that the problems introduced 
by the rendezvous point (e.g. merged or lost tracks) are correctly tackled by all the algorithms. 

Fig. 5.6 shows the OSPA distance of the algorithms. Note again that, in agreement with 
the estimated cardinality distributions, the OSPA distances are nearly identical. 

5.6.2 Scenario 2: 3 TOA sensors 

Figs. 5.7, 5.8 and 5.9 display the statistics (mean and standard deviation) of the estimated 
number of objects obtained, respectively, with M5-GLMB, 5-GLMB and LMB trackers. The 
M5-GLMB and the 5-GLMB tracking filters estimate the object cardinality accurately, while 
the LMB exhibits poor performance and higher standard deviation due to losing some tracks 
when 4 or 5 objects are jointly present in the surveillance area. It is worth noticing that the 
M5-GLMB tracker exhibits performance very close to that of 5-GLMB and that the problems 
introduced by the rendezvous point are again succesfully tackled. 

Fig. 5.10 shows the OSPA distance. Note that the OSPA of M5-GLMB is close to the 
one of 5-GLMB, while LMB shows an overall higher error in agreement with the cardinality 
error due to losing tracks. 
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Time [s] 

Figure 5.6: OSPA distance (c = 600 [m], p = 2) using 1 radar. 
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Figure 5.7: Cardinality statistics for M5-GLMB tracking filter using 3 TOA. 



Figure 5.8: Cardinality statistics for 5-GLMB tracking filter using 3 TOA. 
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The focus of this chapter is on Distributed Multi-Object Tracking (DMOT). We adopt 
the RFS formulation since it provides the concept of probability density for multi-object state 
that allows us to directly extend existing tools in distributed estimation to the multi-object 
case. Such a concept is not available in the Multiple Hypotheses Tracking (MHT) and Joint 
Probabilistic Data Association (JPDA) approaches [Rei79, FS85, FS86, BSF88,BSL95,BP99]. 

DMOF based on the RFS paradigm has previously been described in chapter 4 [Mahl2, 
BCF + 13a, UCJ13, BCF + 14c]. Specifically, an information-theoretic approach to robust dis¬ 
tributed multi-object estimation based on the CCPHD filter has been proposed. However, 
this formulation is not able to provide estimates of the object’s trajectories or tracks. 

Taking into account the above-mentioned considerations, the main contributions in the 
present chapter are [FVV + 15]: 

• the generalization of the robust distributed multi-object filtering approach of subsection 
4.1.1 to multi-object tracking; 

• the development of the Consensus Marginalized 5-Generalized Labeled Multi-Bernoulli 
(CM5GLMB) and Consensus Labeled Multi-Bernoulli (CLMB) filters as the first dis¬ 
tributed multi-object trackers in the RFS framework. 

The proposed solutions are based on the concept of labeled RFS introduced in subsection 
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2.3.6, that enables the estimation of multi-object trajectories in a principled manner [VV13]. 
Moreover, labeled RFS conjugate priors [VV13] have lead to the development of a tractable 
analytic multi-object tracking solution called the 6-Generalized Labeled Multi-Bernoulli (5- 
GLMB) filter [VVP14], However, it is not known if this filter is amenable to DMOT. Nonethe¬ 
less, the M4-GLMB and the LMB filters [FVPV15,RVVD14] are two efficient approximations 
of the d-GLMB filter that 

i. have an appealing mathematical formulation that facilitates an efficient and tractable 
closed-form fusion rule for DMOT; 

ii. preserve key summary statistics of the full multi-object posterior. 

In addition, labeled RFS-based trackers do not suffer from the so-called “ spooky effect ” 
[FSU09] that degrades the performance in the presence of low detection probability like 
the multi-object filters [VVC07,BCF + 13a,UCJ13]. 

6.1 Information fusion with labeled RFS 

In this section, the notion of KLA (4.2) introduced in subsection 4.1.1 will be extended to 
densities of labeled RFSs and will be shown to be equivalent to the normalized weighted 
geometric mean. Proofs of the results are provided in the appendix A. 

6.1.1 Labeled multi-object Kullback-Leibler average 

The benefit of using the RFS framework which provides the concept of probability density of 
the multi-object state, allows us to directly extend the notion of multi-object KLA originally 
devised in [BCF + 13a] to labeled multi-object densities [FVV + 15]. Herewith, we adopt the 
measure theoretic notion of multi-object density given in [VSD05], which does not suffer of 
the unit compatibility problem in the integrals involving the product of powers of multi-object 
densities. Note that this density is equivalent to the FISST density as shown in [VSD05]. 
Thus, the standard inner product notation is extended to multi-object densities 

(f,9) = j /(X)g(X)<SX. (6.1) 

The weighted KLA / KL a °f the labeled multi-object densities f is defined by 

fi<LA - arginf u l D KL (f || /*) , 

* ieN 

J > 0 , w< = !• 

ieN 

where 

Dkl (f li g ) = j /(X) l°g (6-4) 

is the KLD [Mah03, Mah07b] between two multi-object densities /(•) and g(-). 


( 6 . 2 ) 

(6.3) 
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Theorem 3 (KLA of labeled RFSs). 

The weighted KLA defined in (6.2) is the normalized weighted geometric mean of the multi¬ 
object densities /*(•), i G Af, i.e. 

f KLA = © ( w * © /*) • (6-5) 

is AT 

Theorem 3 follows from Theorem 1 of [BCF + 13a, Section III.B]. 

Remark 19. The fusion rule (6.5) is the labeled multi-object version of the multi-object 
Chernoff fusion first proposed by Mahler [MahOO]. 

Subsequently, it will be shown that the KLAs of M5-GLMB and LMB densities are also, 
respectively, a M<5-GLMB and an LMB densities. In particular closed-form solutions for the 
normalized weighted geometric means are derived. These results are necessary to implement 
the fusion rule (6.5) for M5-GLMB and LMB tracking filters. 

6.1.2 Normalized weighted geometric mean of M5-GLMB densities 

The following result holds. 

Theorem 4 (NWGM of M5-GLMB RFSs). 

Let 7r*(-), i = ,T, be M5-GLMB densities on .F(XxIL) and ui l G (0, 1), i = 1 such 

that JffH=i u 1 = 1. Then the normalized weighted geometric mean is given by 

0(u/O tt 1 )(X) = A(X) ^ 5 L (C(N)) [p^]* 

1=1 L&F( L) 

where 

n (#>)"' 

p( L > = _J=3- 

fu (p. (1> )“'dx 
J 1=1 

w< L > = -J—p 

TuLpy’lffiidw.iy'dx] 

FCL i=l L J i=l 

The fusion rule for M<5-GLMBs follows by applying Theorem 4 to find the KLA (6.5) of 
the Md-GLMB’s 7r*(-), i G A f. This is summarized in the following Proposition. 

Proposition 3 (KLA of M<5-GLMB RFSs). 

Suppose that each agent i G Af is provided with an M5-GLMB 7T®(-) and that all agents share 


( 6 . 6 ) 


(6.7) 


( 6 . 8 ) 
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the same label space for the birth process, then the M5-GLMB components p^ L \-) and w ^ 
of the KLA are given by 


where 


P [L \x,i) = 0 (w* ©p (i) )(.T,t) , 
ie N 



i£Af 


(6.9) 

( 6 . 10 ) 

( 6 . 11 ) 

( 6 . 12 ) 


Remark 20. Note that the quantities w ^ and can be independently determined using 
(6.9) and (6.10). Thus, the overall fusion procedure is fully parallelizable. 

Remark 21. Notice that (6.9) is indeed the Cl fusion rule [JU97] for the single-object PDFs. 


6.1.3 Normalized weighted geometric mean of LMB densities 

The following result holds. 


Theorem 5 (NWGM of LMB RFSs). 

Let 7T* = {C 1 = be LMB densities on .F(XxIL) and u l £ (0,1), 

i = 1,... ,T, such that Y%=i N = 1. Then the normalized weighted geometric mean is given 
by 

®K0^) = {(fW,pW)} fgL (6.13) 

1=1 

where 


jJt) 


p(L 


j ft (Ww)"' 

1=1 

n (i-,<«)"'+/n (r;v(x))"dx 

1=1 i=l 


(6.14) 


(6.15) 


The fusion rule for LMBs follows by applying Theorem 5 to find the KLA (6.5) of the 
LMB’s , pP j % £ M. This is summarized in the following Proposition. 

Proposition 4 (KLA of LMB RFSs). 

Suppose that each agent i £ M is provided with an LMB 7r*(-) and that all agents share the 
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same label space for the birth process, then the LMB components j 
are given by 

^ of the KLA 

rW 

_ r 

qW -|- yW 5 

(6.16) 

P {e \x) = 0 (cAOp^Vx) , 

(6.17) 

ieAf 


where 


r^ = j \\ dx , 

i£Af 

(6.18) 

u - ii 0 - AT' • 

(6.19) 

i€Af 



Remark 22. Note that each Bernoulli component can be independently deter¬ 

mined using eqs. (6.16) and (6.17). Thus, the overall fusion procedure is fully parallelizable. 

Remark 23. Notice that also in this case eq. (6.17) is the Cl fusion rule [JU97] for the 
single-object PDFs. 

6.1.4 Distributed Bayesian multi-object tracking via consensus 

At time k, the global KLA (6.5) which would require all the local multi-object densities 
7rj.(-), i G A f, to be available, can be computed in a distributed and scalable way by iterating 
regional averages via the consensus algorithm [BCF + 13a, Section III.A] [BC14] described in 
subsection 2.4.3. Thus, each agent i G M iterates the consensus steps 

*k,l = ® © *1,1- 1) > ( 6 - 2 °) 

jeAf 1 

with 7rj, 0 (-) = TTfc(-); ujhJ © 0) satisfying Jf l jeAf i ul ' 3 = 1> are the consensus weights relating 
agent i to nodes j G A f l . Using the same arguments as the single-object case (subsection 
2.4.3), it follows that, at time k, if the consensus matrix is primitive and doubly stochastic, 
the consensus iterate of each node in the network converges to the global unweighted KLA 
(2.107) of the multi-object posterior densities as the number of consensus iterations l tends 
to infinity. Convergence follows along the same line as in [CA09, BC14] since A"(X x L) is a 
metric space [Mah07b]. In practice, the iteration is stopped at some finite l. 

Remark 24. The consensus iteration (6.20) is the multi-object counterpart of equation 
(2.119) reviewed in Subsection 2.4.3. 

For Md-GLMB multi-object densities, (6.20) can be computed via (6.9) and (6.10), while 
for LMB densities by means of (6.16) and (6.17). In the present work, each single-object 
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density pt is represented by a Gaussian Mixture (GM) of the form 

n g 

p(x ) = ^2 a.j M{x\ Xj, Pj) . (6.21) 

3 = 1 

Note that the fusion rules (6.9) and (6.17) involve exponentiation and multiplication of GM 
in the form of (6.21) which, in general, does not provide a GM. Hence, the same solution 
conceived in subsection 4.1.4 for fusing GMs is also adopted here. 

The other common approach for representing a single object location PDF p is using 
particles. Information fusion involving convex combinations of Dirac delta functions re¬ 
quires, at the best of our knowledge, the exploitation of additional techniques like kernel 
density estimation [UCJ13], least square estimation [HSH + 12, OHD14] or parametric model 
approaches [Coa04] which increase the in-node computational burden. Moreover, the local 
filtering steps are also more resource demanding with respect to a GM implementation. At 
this stage, it is preferred to follow the approach devised in [BCF + 13a] by adopting a GM 
representation. 

6.2 Consensus labeled RFS information fusion 

In this section, two novel fully distributed and scalable multi-object tracking algorithms are 
described by exploiting Propositions 3 and 4 along with consensus [OSFM07, XBL05, CA09, 
BC14] to propagate information throughout the network [FVV + 15]. Pseudo-codes of the 
algorithms are also provided. 

6.2.1 Consensus M5-GLMB filter 

This subsection describes the novel Gaussian Mixture - Consensus Marginalized 5-Generalized 
Labelled Multi-Bernoulli (GM-CMdGLMB) filter algorithm. The operations reported in Table 

6.1 are sequentially carried out locally by each agent i £ M of the network. Each node 
operates autonomously at each sampling interval k, starting from its own previous estimates 
of the multi-object distribution 7r*(-), having location PDFs p^^x,!), W £ I, I £ A"(L), 
represented with a GM, and producing, at the end of the sequence of operations, its new 
multi-object distribution 7T*(-) = Trjy(-) as an outcome of the Consensus procedure. The 
steps of the GM-CdGLMB algorithm follow hereafter. 

1. Each agent i £ Af locally performs a GM-dGLMB prediction (5.31) and update (5.37). 
Additional details of the two procedures can be found in [VVP14, Section IV.B]. 

2. Consensus takes place in each node i involving the subnetwork A/"*. More precisely, 
at each consensus step, node i transmits its data to nodes j such that i £ J\f 3 and 
waits until it receives data from j £ A/" ? \{i}. Next, node i carries out the fusion rule 
of Proposition 3 over A I 1 , i.e. performs (6.20) using local information and information 
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received by Af l . Finally, a merging step for each location PDF is applied to reduce the 
joint communication-computation burden for the next consensus step. This procedure 
is repeatedly applied for a chosen number N > 1 of consensus steps. 

3. After consensus, an estimate of the object set is obtained from the cardinality PMF 
and the location PDFs via an estimate extraction described in Table 6.2. 


Table 6.1: Gaussian Mixture - Consensus Marginalized c)-Generalized Labeled Multi-Bernoulli 
(GM-CM5GLMB) filter 

procedure GM-CM5GLMB(Node i, Time k ) 

Local Prediction > See (5.31) and [VVP14, Table 2, Section V] 

Local Update > See (5.37) and [VVP14, Table 1, Section IV] 

Marginalization > See (5.42) and (5.43) 

for n = 1,..., N do 

Information Exchange 

GM-MdDGLMB Fusion > See eqs. (6.9) and (6.10) 

GM Merging > See [VM06, Table II, Section III.C] 

end for 

Estimate Extraction > See algorithm in Table 6.2 

end procedure 


Table 6.2: GM-CMdGLMB estimate extraction 


Input: 7r 

5 Nmax 


Output: 

X 


for c = 1, 

• • • 5 Cffiax 

do 

P(c ) = 

J2 w 

(I) 


ieM L) 


end for 



C = arg max p(c) 


I = arg max 


I^C Cl 


II 

tj :£e I, 

x = argma xp^(x,l)| 


6.2.2 Consensus LMB filter 

This subsection describes the novel Gaussian Mixture - Consensus Labelled Multi-Bernoulli 
(GM-CLMB) filter algorithm. The operations reported in Table 6.3 will be sequentially 
carried out locally by each agent i 6 Al of the network. Each node operates autonomously 
at each sampling interval k, starting from its own previous estimates of the multi-object 
distribution 7T*(-), having location PDFs p^\ Vf E L, represented with a GM, and producing, 
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at the end of the sequence of operations, its new multi-object distribution 7p(-) = vrj v (-) as an 
outcome of the consensus procedure. A summary description of the steps of the GM-CLMB 
algorithm follows. 

1. Each agent i £ J\f locally performs a GM-LMB prediction (5.59) and update (5.64). The 
update procedure involves two more steps: i) a GM <j-GLMB distribution is created 
from the predicted GM-LMB and updated using the local measurement set Y n) the 
updated GM <j-GLMB is converted back to a GM-LMB distribution. The details of the 
GM-LMB prediction and update can be found in [RVVD14, VVP14], 

2 . Consensus takes place in each node i involving the subnetwork J\f l . More precisely, 
at each consensus step, node i transmits its data to nodes j such that i € M 3 and 
waits until it receives data from j 6 A/"*\{i}. Next, node i carries out the fusion rule 
of Proposition 4 over A I' 1 , i.e. performs (6.20) using local information and information 
received by A f l . Finally, a merging step for each location PDF is applied to reduce the 
joint communication-computation burden for the next consensus step. This procedure 
is repeatedly applied for a chosen number N > 1 of consensus steps. 

3. After consensus, an estimate of the object set is obtained from the cardinality PMF 
and the location PDFs via an estimate extraction described in Table 6.4. 

Table 6.3: Gaussian Mixture - Consensus Labeled Multi-Bernoulli (GM-CLMB) filter 


procedure GM-CLMB(Node i, Time k ) 

Local Prediction > See eq. (5.59) and [RVVD14, Proposition 2, Section III.A] 

GM-LMB -a GM 5-GLMB 
Local Update 

GM 5-GLMB -a- GM-LMB 

> See [RVVD14, Section IV.C.l] 
> See eq. (5.64) and [VVP14, Table 1, Section IV] 
> See [RVVD14, Proposition 4, Section III.B] 

for n = 1 , ..., N do 

Information Exchange 

GM-LMB Fusion 

GM Merging 

end for 

> See eqs. (6.16) and (6.17) 
> See [VM06, Table II, Section III.C] 

Estimate Extraction 

end procedure 

[> See algorithm in Table 6.4 


6.3 Performance evaluation 

To assess performance of the proposed Gaussian Mixture Consensus Marginalized <5-GLMB 
(GM-CM<5GLMB) and LMB (GM-CLMB) described in section 6.2, a 2-dimensional multi¬ 
object tracking scenario is considered over a surveillance area of 50 x 50 [km 2 ], wherein the 
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Table 6.4: GM-CLMB estimate extraction 

Input: 7T = { r(£) >P W } teL > N max 

Output: X 
for c = 1,. .., C max do 
p(c) = w ( L ) 

LeT c ( L) 

end for 

C = arg max p(c) 

L = 0 

for c = 1,..., C do 

L = L U arg max 
£eL\L 

end for 

X = : I G L, x = argmaxp^(x)| 


sensor network of Fig. 6.1 is deployed. The scenario consists of 5 objects as depicted in 
Fig. 6.2. For the sake of comparison, the trackers are also compared with the Gaussian Mix¬ 
ture Consensus CPHD (GM-CCPHD) filter of [BCF + 13a] which, however, does not provide 
labeled tracks. 


xIO 4 



Figure 6.1: Network with 7 sensors: 4 TOA and 3 DOA. 
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xIO 4 



Figure 6.2: Tbject trajectories considered in the simulation experiment. The start/end point 
for each trajectory is denoted, respectively, by «\B. The -k indicates a rendezvous point. 


The kinematic object state is denoted by x = \p x , p x , Py, Py ] T , he. the planar position and 
velocity. The motion of objects is modeled by the filters according to the Nearly-Constant 
Velocity (NCV) model [FS85,FS86,BSF88,BSL95]: 
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( 6 . 22 ) 


where: Wk is a white noise with zero mean and covariance matrix Q, a w = 5 [m/s 2 ] and the 
sampling interval is T s = 5 [s]. 

The sensor network considered in the simulation (see Fig. 6.1) consists of 4 range-only 
(Time Of Arrival, TOA) and 3 bearing-only (Direction Of Arrival, DOA) sensors characterized 
by the following measurement functions: 


{ A(px - x?) + j (p y - y% if i is a DOA sensor 

_ (6.23) 

\J (Px — x 1 ) 2 + ( py — y 1 ) 2 , if i is a TOA sensor 
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where (x 1 , y l ) represents the known position of sensor i. The standard deviation of DOA 
and TOA measurement noises are taken respectively as a do A = 1 [°] and utoa = 100 [to]. 
Because of the non linearity of the aforementioned sensors, the Unscented Kalman Filter 
(UKF) [JU04] is exploited in each sensor in order to update means and covariances of the 
Gaussian components. 

Three different scenarios will be considered, each of which has a different clutter Poisson 
process with parameter A c and a probability of object detection P/> 

• High SNR: A c = 5, Pd = 0.99. These parameters were used in the work [BCF+lSa] 
and, therefore, will be used as a first comparison test. 

• Low SNR: A c = 15, Pd = 0.99. These parameters try to describe, in a realistic way, 
a scenario characterized by high clutter rate A c . 

• Low Pd: A c = 5, Pd = 0.7. These parameters test the distributed algorithms in the 
presence of severe misdetection. 

All the above-mentioned case studies have, for each sensor, a uniform clutter spatial distri¬ 
bution over the surveillance area. 

In the considered scenario, objects pass through the surveillance area with partial prior 
information for object birth locations. Accordingly, a 10-component LMB RFS tvb = 
has been hypothesized for the birth process. Table 6.5 gives a detailed 
summary of such components. Due to the partial prior information on the object birth lo- 

Table 6.5: Components of the LMB RFS birth process at a given time k. 

r W = 0.09 

Pb\x) =Af(x-, mfg , Pb) 

P B = diag(l0 6 ,10 4 ,10 6 ,10 4 ) 


Label 

(k, 1) 

(k, 2 ) 

(k, 3) 

M 

m B 

[0, 0, 40000, 0] T 

[0, 0, 25000, of 

[0, 0, 5000, Of 

Label 

(k, 4) 

(k, 5) 

(k, 6 ) 

(£) 

m B 

[5000, 0, 0, 0] T 

[25000, 0, 0, Of 

[36000, 0, 0, Of 


Label 

(k, 7) 

(k, 8 ) 

(£) 

m B 

[50000, 0, 15000, Of 

[50000, 0, 40000, Of 

Label 

(k, 9) 

(k, 10 ) 

(t) 

m B 

[40000, 0, 50000, Of 

[10000, 0, 50000, of 


cations, some of the LMB components cover a state space region where there is no birth. 
Therefore, clutter measurements are more prone to generate false objects. 

Multi-object tracking performance is evaluated in terms of the Optimal SubPattern As¬ 
signment (OSPA) metric [SVV08] with Euclidean distance p = 2 and cutoff c = 600 [to]. The 
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reported metric is averaged over 100 Monte Carlo trials for the same object trajectories but 
different, independently generated, clutter and measurement noise realizations. The duration 
of each simulation trial is fixed to 1000 [s] (200 samples). 

The GM-CMdGLMB and the GM-CLMB are capped to 20000, 8000 and 3000 hypotheses 
[VV13, VVP14], respectively, for the High SNR, Low SNR and Low Pd scenario, and are 
coupled with the parallel CPHD look ahead strategy described in [VV13, VVP14]. The CPHD 
filter is similarly capped, for each case study, to the same number of components through 
pruning and merging of mixture components. 

The parameters of the GM-CCPHD filter have been chosen as follows: the survival prob¬ 
ability is P s = 0.99; the maximum number of Gaussian components is N max = 25; the 
merging threshold is = 4; the truncation threshold is 7 1 = 10 , the extraction threshold 

is 7 e = 0.5; the birth intensity function is the PHD of the LMB RFS of Table 6.5. A single 
consensus step L = 1 is employed for all the simulations. 

6.3.1 High SNR 

Figs. 6.3, 6.4 and 6.5 display the statistics (mean and standard deviation) of the esti¬ 
mated number of objects obtained, respectively, with GM-CCPHD, GM-CLMB and GM- 
CMdGLMB. As it can be seen, all the distributed algorithms estimate the object cardinality 
accurately, with the GM-CM5GLMB exhibiting better estimated cardinality variance. Note 
that the difficulties introduced by the rendezvous point (e.g. merged or lost tracks) are 
correctly tackled by all the distributed algorithms. 

Fig. 6.6 shows the OSPA distance for the three algorithms. The improved localization 
performance of GM-CLMB and GM-CMdGLMB is attributed to two factors: (a) the “spooky 
effect” [FSU09] causes GM-CCPHD filter to temporarily drop objects which are subjected to 
missed detections and to declare multiple estimates for existing tracks in place of the dropped 
objects, and (b) the two trackers are generally able to better localize objects due to a more 
accurate propagation of the posterior density. Note that GM-CLMB and GM-CMdGLMB 
have similar performance since the approximations introduced by the LMB tracker (see (5.65)- 
(5.66)) are negligible in the case of high SNR. 

6.3.2 Low SNR 

Figs. 6.7 and 6.8 display the statistics (mean and standard deviation) of the estimated 
number of objects obtained, respectively, with the GM-CCPHD and the GM-CM5GLMB. 
As it can be seen, the two distributed algorithms estimate the object cardinality accurately, 
with the GM-CMdGLMB exhibiting again better estimated cardinality variance. 

Note that the GM-CLMB fails to track the objects. The problem is due to the approxi¬ 
mation (5.65)-(5.66) made to convert a d-GLMB to an LMB, becoming significant with low 
SNR. In particular, each local tracker fails to properly set the existence probability of the 
tracks for three main factors: (a) no local observability, (b) high clutter rate and (c) loss of 
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Figure 6.3: Cardinality statistics for the GM-CCPHD filter under high SNR. 
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the full posterior cardinality distribution after the probability density conversion. By having 
low existence probabilities, the extraction of the tracks fails even if the single object densities 
are correctly propagated in time. 

Fig. 6.9 shows the OSPA distance for the current scenario. As for the previous case study, 
GM-CM5GLMB outperforms GM-CCPHD. 

6.3.3 Low Pu 

Figs. 6.10 displays the statistics (mean and standard deviation) of the estimated number of 
objects obtained with the GM-CMdGLMB. It is worth pointing out that the only working 
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Figure 6.5: Cardinality statistics for the GM-CMAGLMB tracker under high SNR. 



Figure 6.6: OSPA distance (c = 600 [m], p = 2) under high SNR. 


distributed algorithm is, indeed, the GM-CM(5GLMB and that it exhibits good performance 
in terms of average number of estimated objects in a very tough scenario with Pd = 0.7. 
Fig. 6.11 shows the OSPA distance for the current scenario. 
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Figure 6.7: Cardinality statistics for the GM-CCPHD filter under low SNR. 
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^ure 6.8: Cardinality statistics for the GM-CM5GLMB tracker under low SNR. 



Figure 6.9: OSPA distance (c = 600 [m], p = 2) under low SNR. 
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Figure 6.11: OSPA distance (c = 600 [m], p = 2) under low Pp. 
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Conclusions and future 

work 


The present dissertation has summarized the work carried out during my 3-year Ph.D. re¬ 
search work in the field of distributed filtering and tracking. Contributions have spanned 
from single-object to multi-object scenarios, in both centralized and distributed settings. In 
particular, distributed tracking over sensor networks requires redesigning the architecture 
and algorithms to address the following issues: 

1. lack of a central fusion node; 

2. scalable processing with respect to the network size; 

3. each node operates without knowledge of the network network topology; 

4. each node operates without knowledge of the dependence between its own information 
and the information received from other nodes. 

main ingredients have been considered in order to provide novel research contributions. 

The Bayesian framework for modeling the quantities to be estimated as random vari¬ 
ables characterized by Probability Density Functions (PDFs), and for providing an 
improved estimation of such quantities by conditioning the PDF on the available noisy 
measurements. 

ii. The Random Finite Set (RFS) formulation as it provides the concept of probability 
density of the multi-object state that allows us to directly generalize single-object esti¬ 
mation to the multi-object case. 

hi. The Kullback-Leibler Average (KLA) as a sensible information-theoretic definition for 
consistently fusing PDFs of the quantity to be estimated provided by different nodes 
(or agents). 

iv. Consensus as a tool for developing scalable and reliable distributed fusion techniques 
over a network where each agent aims to compute the collective average of a given 
quantity by iterative regional averages, where the terms “collective” and “regional” 
mean “over all network nodes” and, respectively, “over neighboring nodes only”. 


Four 

i. 
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The main contributions of the present thesis are herewith summarized. 

□ A novel distributed single-object filter, namely Parallel Consensus on Likelihoods and 
Priors (CLCP). The proposed algorithm is based on the idea of carrying out, in par¬ 
allel, a separate consensus for the novel information (likelihoods) and one for the prior 
information (priors). This parallel procedure is conceived as an improvement of the 
Consensus on Posteriors (CP) approach to avoid underweighting the novel information 
during the fusion steps. The outcomes of the two consensuses are then combined to 
provide the fused posterior density. 

□ Two novel consensus Multiple Model (MM) filters to be used for tracking a maneuvering 
object, namely Distributed First Order Generalized Pseudo-Bayesian (DGPBi) and 
Distributed Interacting Multiple Model (DIMM) filters. 

□ A novel result on Distributed Multi-Object Filtering (DMOF) over a sensor network 
using RFSs. Specifically, a generalisation of the single-object Kullback-Leibler Average 
(KLA) to the RFS framework has been devised. 

□ A novel consensus multi-object filter, the Consensus Cardinalized Probability Hypothesis 
Density (CCPHD) filter. Each tracking agent locally updates multi-object CPHD, i.e. 
the cardinality distribution and the PHD, exploiting the multi-object dynamics and the 
available local measurements, exchanges such information with communicating agents 
and then carries out a fusion step to combine the information from all neighboring 
agents. 

□ A novel approximation of the d-GLMB filter, namely the Marginalized 5-Generalised La¬ 
beled Multi-Bernoulli (Md-GLMB). The result is based on a principled Generalised La¬ 
beled Multi-Bernoulli (GLMB) approximation of the labeled RFS posterior that matches 
exactly the posterior PHD and cardinality distribution. The proposed approximation 
can be interpreted as performing a marginalization with respect to the association his¬ 
tories arising from the 5-GLMB filter. The key advantage of the new filter lies in the 
reduced growth rate of the number of new components generated at each filtering step. 
In particular, the approximation (or marginalization) step performed after each update 
is guaranteed to reduce the number of generated components which normally arise from 
multiple measurement-to-track association maps. Typically, the proposed M5-GLMB 
filter requires much less computation and storage especially in multi-sensor scenarios 
compared to the 5-GLMB filter. Furthermore the proposed Md-GLMB filter inherits 
the same implementation strategies and parallelizability of the d-GLMB filter. 

□ The generalisation of the KLA to labeled RFSs. 

□ Two novel consensus tracking filters, namely a Consensus Marginalized 5-Generalized 
Labeled Multi-Bernoulli (CMJGLMB) and Consensus Labeled Multi-Bernoulli (CLMB) 
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tracking filter. The proposed algorithms provide fully distributed, scalable and compu¬ 
tationally efficient solutions for multi-object tracking. 

From the above-mentioned contributions, the following conclusions can be drawn. 

> The KLA turns out to be an effective information-theoretic approach for distributing 
information over a sensor network. 

> Consensus algorithms can be exploited so as to fuse, in a fully distributed and scal¬ 
able fashion, the information collected from multiple sensors which are supposed to be 
heterogeneous and geographically dispersed. 

> The CLCP has proven to be a reliable distributed filter. Its performance is crucially de¬ 
pendant on the choice of a suitable weight of the innovation (likelihood) term combined 
with the prior. In particular, the choice of such a weight equal to the (usually unknown) 
number of nodes provides the best performance when the number of consensus steps is 
sufficiently large, but for very few consensus steps can be improved resorting to other 
simple choices, proposed in this thesis work, which do not require knowledge of the 
network size. 

> The DGPBi and DIMM have been proven to be capable of tracking a highly-maneuvering 
object with a network consisting of range-only and angle-only sensors along with com¬ 
munication nodes. DIMM performs significantly better than DGPBi when a turn is 
approaching. 

> The multi-object KLA in the RFS framework turns out to be a significant contribution 
in the field of distributed multi-object filtering. The proposed consensus algorithm 
admits an intuitive interpretation in terms of distributed averaging of the local multi¬ 
object densities. Its formulation has led to the development of the CCPHD filter which 
proved its effectiveness in a realistic simulation scenario. 

> The labeled RFS approach turns out to be a strong mathematical tool for modeling the 
probability density of a multi-object state. Its key features allow to directly generalise, 
in a rigorous and elegant way, the multi-object KLA to the broader labeled multi-object 
case. 

> The Md-GLMB provides an excellent tradeoff between computational/storage require¬ 
ments and performance for multi-object tracking. Moreover, it has an appealing math¬ 
ematical formulation which 1) facilitates an efficient and tractable closed-form solution 
for distributed multi-object tracking and 2) preserves key summary statistics of the full 
multi-object posterior. 

i> The CLMB fails to track objects in low SNR scenarios. The problem is due to the 
approximation made to convert a d-GLMB into an LMB, becoming significant with 
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CHAPTER 7. CONCLUSIONS AND FUTURE WORK 


low SNR. In particular, each local tracking filter fails to properly set the existence 
probability of the tracks in the case of having: 1) no local observability, 2) high clutter 
rate and 3) loss of the full posterior cardinality distribution after the probability density 
conversion. By having low existence probabilities, the extraction of the tracks fails even 
if the single object densities are correctly propagated in time. 

Based on the above achievements, possible future work will concern the following issues. 

o To consider sensors with different field-of-view to cope with real-world sensor capabili¬ 
ties. So far, it has been assumed that the area subjected to surveillance is exactly the 
intersection of the sensor fields-of-view, or that each sensor field-of-view is exactly the 
surveillance region. By a practical viewpoint this issue needs to be addressed. Thus, 
the proposed algorithms will have to take into account a state-dependent probability 
of detection instead of a constant one as assumed in the present work. Moreover, the 
fusion techniques will have to be redesigned to consider the above-mentioned context. 

o To investigate distributed measurement-driven object initialization. The object birth 
process of the RFS filtering approach is capable of accounting for a large variety of sce¬ 
narios in which object birth locations could also be unknown. However, a measurement- 
driven approach would be much more effective compared to a static birth process. The 
information provided by the sensors could be used to dynamically accommodate each 
sensor birth process. Thus, a more complete picture of newborn objects could be 
achieved by fusing local information with the one provided by the other nodes. 
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Proofs 



Proof of Theorem 1: 

The cost to be minimized by the KLA is 


J (m) = E nl Dkl (h p) = E ^ E ^ lo & 
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where the relationship Y^iejV tt* = 1 has been exploited. Let /j be defined as in (3.24), then 


J (h) = E h j iog( j=j 


jen 


c 


(A.l) 


where c — E II (Py • Hence 
hen ieM 


J (n) = D i<l (/i||//)- log c 


(A.2) 


where the relationship YljeTZ /H = 1 has been exploited. Then the infimum of J (/t), i.e. the 
weighted KLA of (3.23), is just provided by // = ~p as in (3.24). □ 
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Proof of Theorem 2: 

Let 


w) = n [/*« 

i 


Wi 


C = 



(A.3) 

(A.4) 


so that the fused multi-object density in (4.4) can be expressed as f(X ) = f(X)/c. Then, 
the cost to be minimized by the KLA is 


j(f) = YjhDkl (/ I! f) 

i 




= D kl (/ || /) - log(c) 


Since the KLD is always nonnegative and is zero if and only if its two arguments coincide 
almost everywhere, the above J(f) is trivially minimized by f kla(X) defined as in (4.4). □ 


Proof of Proposition 1: 

First, the cardinality distributions of 7r(-) and tt(-) are proven to be the same. Following 
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(2.81), the cardinality distribution of the GLMB in (5.23) is given by: 


Pin) = Y w [I \L) (A.5) 

LeT n ( L) l€T{h) 

= Y w ( L ) ( A - 6 ) 

Ler n ( l) 

=-w Z! (A- 7 ) 

n ' (H,...bt)eL« 

= — Y Tr({(xi,£i),...,(x n ,£ n )})dxi---dx n (A.8) 

n ' (4,...G)6L«' / 

= — / 7r({xi,..,.,x n )})dx 1 ---dx ri (A.9) 

n! 7 

= p(n) (A.10) 

Second, the PHD corresponding to 7r(-) and 7f(-) are proven to be the same. Note from the 
definition of p^\x, l) that 

p ({e,e i-.M)( X j £) = 1 {£A) ... A} (^)p {4; a} (x,£) (A.11) 

= J p({(x,£),(xi,£i),...,(x n ,£ n )})dxi---dx n (A.12) 


Substituting the GLMB density (5.23) into the PHD of the object with label (or track) £ 


yields 

/» 


d(x,£) = / 7r({(x, £) U X}) <5X 

from which 

(A.13) 
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(A.15) 

= Y Y ({£} U L) p^\-,£)dx p( J \x,£) 

ie j’(L) ie^(L) 

(A.16) 

where the last step follows from Lemma 3 in [VV13, Section III.B], Noting that p^\-,£) is 
a probability density, and using (5.24) gives 

d(x,£)= Y Y 5 I ({£}UL)w{I)p (I \x,£) 

LeT( L) IeT(L) 

(A.17) 

= Yj t»({f}Ui)p*^(i,l) 
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from which by applying (A.12), we obtain: 


A 00 1 _ r 

d ( x ^) w ^ w {{£Ji ,. • . 4 }) J p({{x,£), (xi,£i),{x n ,£ n )}) dxi ■ ■ ■ dx r 


n=0 ■ (7i,...,7 n )eL' 


= J tt({(x,.£)}UX)AX 

= d(x, £). 


(A.20) 
(A.21) 
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□ 


Proof of Proposition 2: 

We apply the result in Proposition 1 which can be used to calculate the parameters of the 
marginalized 4GLMB density. Notice that such a result applies to any labeled RFS density 
and our first step is to rewrite the 5-GLMB density (2.86) in the general form for a labeled 
RFS density [PVV+14], i.e. 

7r (X) = w(£(X))p(X) (A.23) 

where 
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(A.29) 


Applying Proposition 1, the parameters and p^\-) for the M4GLMB approximation 
that match the cardinality and PHD are 


w {i \l) = 4(A)w(J) = 4(4 E w{1)0 (A.30) 

«eE 
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and 


P (I \x ,£) = 1 i{£)pi-{e}{x,£) (A.31) 

= 1/(7) Jp({(x,£),(xi,£i),...,(x j ,£j)})d(xi,...,x j ) (A.32) 

where I\{£} = {£\,... ,£j}- Substituting the expression (A.29) in (A.32), we have 
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«GH 


and, noting that / = {7, £\,... ,£j}, it follows that only one term in the sum over J is non-zero 
thus giving 

P {I) (x,£) = lj(l)A({(x,l),(xi,li (xj,£j)}) 1 E w{I '°P ( ^\x,£) (A.35) 

E W €GH 
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Consequently, the IvM-GLMB approximation is given by 
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Proof of Theorem f: 

For the sake of simplicity, let us consider only two Md-GLMB densities i £ {1, 2}. From (6.6) 
one gets 



(A.42) 


Notice that the exponentiation of a sum of delta functions is a sum of the exponentiated 
delta function terms, i.e 

A(X) (y Ll) ) W [(p (il) ) W ] X if £(X) = U 

7f(xr = < : : 

A(X) [( p^Y ] X if £(X) = L n 

Thus, (A.42) yields: 
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L 
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The normalization constant K can be easily evaluated exploiting Lemma 3 of [VV13, Section 
III.B], i.e. 
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Applying (A.46) in (A.44) one has 
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It can be proved by induction that Theorem 4 holds considering I Md-GLMB densities instead 
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of 2. 


□ 


Proof of Proposition 3: 

The proof readily follows by noting that the KLA (6.5) can be evaluated using (6.6) of 
Theorem 4. □ 


The following Lemma is useful to prove Theorem 5. 

Lemma 1 (Normalization constant of LMB KLA). 

Let 7r 1 (X) = | | and 7t 2 (X) = } be ^ wo LMB densities on X x L 

and uj £ (0,1) then 


K = j 7r 1 (X) tJ 7r 2 (X) 1 ~ a; 5X 
= (uj © 7T 1 , (1 - Uj) © 7 r 2 j 

= (V') + . 

Proof of Lemma 1: 

We make use of the Binomial Theorem [AS64] which states 

E/ i = (i + /) L 

LCL 

Applying Lemma 3 of [VV13, Section III.B] to the definition (A.49) of K gives 
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Applying (A.52) to (A.53) gives 
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□ 


Proof of Theorem 5: 

For the sake of simplicity, let us consider only two LMB densities i £ {1,2}. From (6.13) one 
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gets 
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Hence, recalling definitions (A.56) and (A.57), one has 
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It can be proved by induction that Theorem 5 holds considering I LMB densities instead of 

2 . □ 


Proof of Proposition f: 

The proof readily follows by noting that the KLA (6.5) can be evaluated using (6.13) of 
Theorem 5. □ 


134 



Bibliography 


[AF70] 

[Aka98] 

[AM12] 

[AS64] 

[BBS88] 

[BC14] 

[BCF+13a] 

[BCF+13b] 

[BCF+13c] 

[BCF+13d] 

[BCF14a] 


G. Ackerson and K. S. Fu, “On state estimation in switching environments,” 
IEEE Trans, on Automatic Control, vol. 15, no. 1, pp. 10-17, 1970. 

H. Akaike, “Information theory and the extension of the maximum likelihood 
principle,” in Selected Papers of Hirotugu Akaike. Springer New York, 1998, pp. 
199-213. 

B. Anderson and J. B. Moore, Optimal filtering. Courier Dover Publications, 

2012. 

M. Abramowitz and I. A. Stegun, Handbook of mathematical functions: with 
formulas, graphs, and mathematical tables. Washington, DC: National Bureau 
of Standards, 1964. [Online]. Available: http://people.math.sfu.ca/~cbm/aands/ 
abr aniowi t z _and _s t egun. p df 

H. Blorn and Y. Bar-Shalom, “The interacting multiple model algorithm for sys¬ 
tems with Markovian switching coefficients,” IEEE Trans, on Automatic Control, 
vol. 33, no. 8, pp. 780-783, 1988. 

G. Battistelli and L. Chisci, “Kullback-Leibler average, consensus on probability 
densities, and distributed state estimation with guaranteed stability,” Automat- 
ica, vol. 50, no. 3, pp. 707-718, 2014. 

G. Battistelli, L. Chisci, C. Fantacci, A. Farina, and A. Graziano, “Consensus 
CPHD filter for distributed multitarget tracking,” IEEE Journal of Selected Top¬ 
ics in Signal Processing, vol. 7, no. 3, pp. 508-520, 2013. 

—, “Distributed multitarget tracking with range-Doppler sensors,” in Proc. 
21 st European Signal Processing Conference (EUSIPCO), pp. 1-5, Marrakech, 
Morocco, 2013. 

—, “Distributed multitarget tracking for passive multireceiver radar systems,” 
in Proc. lfi h International Radar Symposium (IRS), pp. 337-342, Dresden, Ger¬ 
many, 2013. 

—, “Distributed tracking with Doppler sensors,” in Proc. 52 nd IEEE Conf. on 
Decision and Control (CDC), pp. 4760-4765, Florence, Italy, 2013. 

G. Battistelli, L. Chisci, and C. Fantacci, “Parallel consensus on likelihoods and 
priors for networked nonlinear filtering,” IEEE Signal Processing Letters, vol. 21, 
no. 7, pp. 787-791, 2014. 


135 



BIBLIOGRAPHY 


[BCF+14b] 

[BCF+14c] 

[BCMP11] 

[BJA12] 

[BP99] 

[BSF88] 

[BSL95] 

[BSLK01] 

[BVV14] 

[CA09] 

[Cam70] 

[CCMIO] 

[CCSZ08] 

[CJMRIO] 


G. Battistelli, L. Chisci, C. Fantacci, A. Farina, and A. Graziano, “Consensus- 
based multiple-model Bayesian filtering for distributed tracking,” IET Radar, 
Sonar & Navigation, DOI: 10.1049/iet-rsn.2014.0071, 2014. 

G. Battistelli, L. Chisci, C. Fantacci, N. Forti, A. Farina, and A. Graziano, “Dis¬ 
tributed peer-to-peer multitarget tracking with association-based track fusion,” 
in Proc. 17 th Int. Conf. on Information Fusion (FUSION), Salamanca, Spain, 
2014. 

G. Battistelli, L. Chisci, S. Morrocchi, and F. Papi, “An information-theoretic 
approach to distributed state estimation,” in Proc. 18 th IFAC World Congress, 
pp. 12477-12482, Milan, Italy, 2011. 

T. Bailey, S. Julier, and G. Agamennoni, “On conservative fusion of information 
with unknown non-Gaussian dependence,” in Proc. 15 th Int. Conf. on Informa¬ 
tion Fusion (FUSION), pp. 1876-1883, Singapore , 2012. 

S. Blackman and R. Popoli, Design and analysis of modern tracking systems. 
Artech House, Norwood, MA, USA, 1999. 

Y. Bar-Shalom and T. Fortmann, Tracking and data association. Academic 
Press, San Diego, CA, USA, 1988. 

Y. Bar-Shalom and X. Li, Multitarget-multisensor tracking: principles and tech¬ 
niques. YBS Publishing, Storrs, CT, USA, 1995. 

Y. Bar-Shalom, X. R. Li, and T. Kirubarajan, Estimation with applications to 
tracking and navigation. John Wiley & Sons, Hoboken, NJ, USA, 2001. 

M. Beard, B.-T. Vo, and B.-N. Vo, “Multi-target tracking with merged measure¬ 
ments using labelled random finite sets,” in Proc. 17 th Int. Conf. on Information 
Fusion (FUSION), Salamanca, Spain, 2014. 

G. C. Calafiore and F. Abrate, “Distributed linear estimation over sensor net¬ 
works,” International Journal of Control, vol. 82, no. 5, pp. 868-882, 2009. 

L. Campbell, “Equivalence of Gauss’s principle and minimum discrimination 
information estimation of probabilities,” The Annals of Mathematical Statistics, 
vol. 41, no. 3, pp. 1011-1015, 1970. 

K. C. Chang, C.-Y. Chong, and S. Mori, “Analytical and computational evalu¬ 
ation of scalable distributed fusion algorithms,” IEEE Trans, on Aerospace and 
Electronic Systems, vol. 46, no. 4, pp. 2022-2034, 2010. 

R. Carli, A. Chiuso, L. Schenato, and S. Zampieri, “Distributed Kalman filtering 
based on consensus strategies,” IEEE Journal on Selected Areas in Communica¬ 
tions, vol. 26, pp. 622-633, 2008. 

D. Clark, S. Julier, R. Mahler, and B. Ristic, “Robust multi-object fusion with 
unknown correlation,” in Proc. 2010 Sensor Signal Processing for Defense Con¬ 
ference (SSPD 2010), London, UK, 2010. 


136 



BIBLIOGRAPHY 


[CMC90] 

[Coa04] 

[CS10] 

[CT12] 

[FBC+12] 

[FFTS10] 

[Fisl2] 

[FS85] 

[FS86] 

[FSU09] 

[FVPV15] 

[FVV+15] 

[Gau04] 

[GMN97] 

[GSW09] 


C.-Y. Chong, S. Mori, and K.-C. Chang, “Distributed multitarget multisensor 
tracking,” in Multitarget-Multisensor Tracking: Advanced Applications , Y. Bar- 
Shalorn, Ed. Artech House, Norwood, MA, USA, 1990, ch. 8. 

M. Coates, “Distributed particle filters for sensor networks,” in Proc. 3 rd Inter¬ 
national Symposium on Information Processing in Sensor Networks (IPSN'Of), 
pp. 99-107, New York, NY, USA, 2004. 

F. S. Cattivelli and A. Sayed, “Diffusion strategies for distributed Kalman filter¬ 
ing and smoothing,” IEEE Trans, on Automatic Control, vol. 55, pp. 2069-2084, 
2010. 

T. M. Cover and J. A. Thomas, Elements of information theory. John Wiley & 
Sons, Hoboken, NJ, USA, 2012. 

C. Fantacci, G. Battistelli, L. Chisci, A. Farina, and A. Graziano, “Multiple- 
model algorithms for distributed tracking of a maneuvering target,” in Proc. 15 th 
Int. Conf. on Information Fusion (FUSION), Singapore, 2012, pp. 1028-1035. 

M. Farina, G. Ferrari- 1 Trecate, and R. Scattolini, “Distributed moving horizon 
estimation for linear constrained systems,” IEEE Trans, on Automatic Control, 
vol. 55, no. 11, pp. 2462-2475, 2010. 

R. A. Fisher, “On an absolute criterion for fitting frequency curves,” Messenger 
of Math., pp. 41-155, 1912. 

A. Farina and F. A. Studer, Radar data processing, vol. I: introduction and 
tracking. Research Studies Press, Letchworth, Hertfordshire, England, 1985. 

—, Radar data processing, vol. II: advanced topics and applications. Research 
Studies Press, Letchworth, Hertfordshire, England, 1986. 

D. Franken, M. Schmidt, and M. Ulrnke, ““Spooky action at a distance” in the 
cardinalized probability hypothesis density filter,” IEEE Trans, on Aerospace 
and Electronic Systems, vol. 45, no. 4, pp. 1657-1664, 2009. 

C. Fantacci, B.-T. Vo, F. Papi, and B.-N. Vo, “The marginalized J-GLMB filter,” 
preprint available online at arXiv:1501.00926, 2015. 

C. Fantacci, B.-N. Vo, B.-T. Vo, G. Battistelli, and L. Chisci, “Consensus la¬ 
beled random finite set filtering for distributed multi-object tracking,” preprint 
available online at arXiv:1501.01579, 2015. 

F. C. Gauss, Theoria Motus Corporum Coelestium in Sectionibus Conicis Solem 
Ambientium. Hamburg, 1809. (Translation of Theoria Motus: Theory of the 
Motion of the Heavenly Bodies Moving About the Sun in Conic Sections. Dover 
Phoenix ed., 2004). 

I. Goodman, R. P. S. Mahler, and H. Nguyen, Mathematics of Data Fusion. 
Kluwer Academic Publishers, Norwell, MA, USA, 1997. 

R. Georgescu, S. Schoenecker, and P. Willett, “GM-CPHD and MLPDA applied 
to the SEABAR07 and TNO-blind multi-static sonar data,” in Proc. 12 th Int. 
Conf. on Information Fusion (FUSION), pp. 1851-1858, Seattle, WA, USA, 2009. 


137 



BIBLIOGRAPHY 


[HL64] 

[HSH+12] 

[Hur02] 

[Jay03] 

[Jaz07] 

[JBU06] 

[JU97] 

[JU04] 

[JUDW95] 

[Jul06] 

[Jul08] 

[Kal60] 

[KB61] 

[Kol50] 

[KT08] 


Y. Ho and R. Lee, “A Bayesian approach to problems in stochastic estimation 
and control,” IEEE Trans, on Automatic Control , vol. 9, no. 4, pp. 333-339, 
1964. 

O. Hlinka, O. Sluciak, F. Hlawatsch, P. M. Djuric, and M. Rapp, “Likelihood 
consensus and its application to particle filtering,” IEEE Trans, on Signal Pro¬ 
cessing , vol. 60, no. 8, pp. 4334-4349, 2012. 

M. Hurley, “An information theoretic justification for covariance intersection and 
its generalization,” in Proc. 5 th Int. Conf. on Information Fusion (FUSION), pp. 
505-511, Annapolis, MD, USA, 2002. 

E. T. Jaynes, Probability Theory: The Logic of Science. Cambridge University 
Press, 2003. 

A. H. Jazwinski, Stochastic processes and filtering theory. Courier Dover Pub¬ 
lications, 2007. 

S. J. Julier, T. Bailey, and J. K. Uhlmann, “Using exponential mixture models for 
suboptimal distributed data fusion,” in Proc. 2006 IEEE Nonlinear Stat. Signal 
Proc. Workshop (NSSPW'06), pp. 160-163, Cambridge, UK, 2006. 

S. J. Julier and J. K. Uhlmann, “A non-divergent estimation algorithm in the 
presence of unknown correlations,” in Proc. IEEE American Control Conference 
(ACC 1997), pp. 2369-2373, Albuquerque, NM, 1997. 

—, “Unscented filtering and nonlinear estimation,” Proc. of the IEEE, vol. 92, 
no. 3, pp. 401-422, 2004. 

S. J. Julier, J. K. Uhlmann, and H. F. Durrant-Whyte, “A new approach for 
filtering nonlinear systems,” in Proc. 1995 American Control Conference (ACC), 
1995, pp. 1628-1632, Seattle, WA, USA, 1995. 

S. J. Julier, “An empirical study into the use of Chernoff information for ro¬ 
bust, distributed fusion of Gaussian mixture models,” in Proc. 9 th Int. Conf. on 
Information Fusion (FUSION), Florence, Italy, 2006. 

—, “Fusion without independence,” 2008 IET Seminar on Target Tracking 
and Data Fusion: Algorithms and Applications, 2008. 

R. E. Kalman, “A new approach to linear filtering and prediction problems,” 
Trans. ASME - Journal of Basic Engineering, vol. 82 (Series D), no. 1, pp. 35- 
45, 1960. 

R. E. Kalman and R. S. Bucy, “New results in linear filtering and prediction 
theory,” Trans. ASME - Journal of Basic Engineering, vol. 83 (Series D), no. 1, 
pp. 95-108, 1961. 

A. N. Kolmogorov, Foundations of the Theory of Probability. New York, Chelsea 
Pub. Co., 1950. 

M. Kamgarpour and C. Tomlin, “Convergence properties of a decentralized 
Kalman filter,” in Proc. 47 th IEEE Conf. on Decision and Control (CDC), pp. 
3205-3210, Cancun, Mexico, 2008. 


138 



BIBLIOGRAPHY 


[Lee08] 

[LJ05] 

[LJ12] 

[Mah36] 

[MahOO] 

[Mah03] 

[Mah04] 

[Mah07a] 

[Mah07b] 

[Mahl2] 

[Mahl3] 
[Mahl4] 
[Mat 75] 
[MSS62] 

[OCIO] 

[OHD14] 


D. J. Lee, “Nonlinear estimation and multiple sensor fusion using unscented 
information filtering,” IEEE Signal Processing Letters, vol. 15, no. 8, pp. 861- 
864, 2008. 

X. R. Li and V. Jilkov, “Survey of maneuvering target tracking. Part V: Multiple- 
model methods,” IEEE Trans, on Aerospace and Electronic Systems, vol. 41, 
no. 4, pp. 1255-1321, 2005. 

W. Li and Y. Jia, “Consensus-based distributed multiple model UKF for jump 
Markov nonlinear systems,” IEEE Trans, on Automatic Control, vol. 57, no. 1, 
pp. 227-233, 2012. 

P. C. Mahalanobis, “On the generalised distance in statistics,” Proc. National 
Institute of Sciences of India, vol. 2, no. 1, pp. 49-55, 1936. 

R. P. S. Mahler, “Optimal/robust distributed data fusion: a unified approach,” 
Proc. SPIE Defense and Security Symposium 2000, 2000. 

—, “Multitarget Bayes filtering via first-order multitarget moments,” IEEE 
Trans, on Aerospace and Electronic Systems, vol. 39, no. 4, pp. 1152-1178, 2003. 

—, ““Statistics 101” for multisensor, multitarget data fusion,” IEEE Trans, 
on Aerospace and Electronic Systems, vol. 19, no. 1, pp. 53-64, 2004. 

—, “PHD filters of higher order in target number,” IEEE Trans, on Aerospace 
and Electronic Systems, vol. 43, no. 4, pp. 1523-1543, 2007. 

—, Statistical multisource-multitarget information fusion. Artech House, Nor¬ 
wood, MA, 2007. 

—, “Toward a theoretical foundation for distributed fusion,” in Distributed 
Data Fusion for Network-Centric Operations, D. Hall, M. L. II, C.-Y. Chong, 
and J. Llinas, Eds. CRC Press, Boca Raton, FL, USA, 2012, ch. 8. 

—, ““Statistics 102” for multisource-multitarget detection and tracking,” IEEE 
Journal of Selected Topics in Signal Processing, vol. 7, no. 3, pp. 376-389, 2013. 

—, Advances in Statistical Multisource-Multitarget Information Fusion. 
Artech House, Norwood, MA, 2014. 

G. Matheron, Random sets and integral geometry. John Wiley & Sons, Hoboken, 
NJ, USA, 1975. 

L. A. McGee, S. F. Schmidt, and G. L. Smith, “Applications of statistical filter 
theory to the optimal estimation of position and velocity on board a circumlunar 
vehicle,” NASA Technical Report R-135, Tech. Rep., 1962. 

B. N. Oreshkin and M. J. Coates, “Asynchronous distributed particle filter via 
decentralized evaluation of Gaussian products,” in Proc. 13 th Int. Conf. on In¬ 
formation Fusion (FUSION), Edinburgh, UK, 2010. 

F. H. O. Hlinka and P. M. Djuric, “Consensus-based distributed particle filtering 
with distributed proposal adaptation,” IEEE Trans, on Signal Processing, vol. 62, 
no. 12, pp. 3029-3041, 2014. 


139 



BIBLIOGRAPHY 


[OS07] R. Olfati-Saber, “Distributed Kalman filtering for sensor networks,” in Proc. 46 th 
IEEE Conf. on Decision and Control (CDC'), pp. 5492-5498, New Orleans, LA, 
USA, 2007. 

[OSFM07] R. Olfati-Saber, J. Fax, and R. Murray, “Consensus and cooperation in networked 
multi-agent systems,” Proceedings of the IEEE , vol. 95, no. 1, pp. 49-54, 2007. 

[PVV + 14] F. Papi, B.-N. Vo, B.-T. Vo, C. Fantacci, and M. Beard, “Generalized labeled 
multi-Bernoulli approximation of multi-object densities,” preprint available on¬ 
line at arXiv.1412.5294-, 2014. 

[RAG04] B. Ristic, S. Arulampalm, and N. J. Gordon., Beyond the Kalman filter: Particle 
filters for tracking applications. Artech House Publishers, 2004. 

[Rei79] D. Reid, “An algorithm for tracking multiple targets,” IEEE Trans, on Automatic 
Control , vol. 24, no. 6, pp. 843-854, 1979. 

[RVVD14] S. Reuter, B.-T. Vo, B.-N. Vo, and K. Dietrmayer, “The labeled multi-Bernoulli 
filter,” IEEE Trans, on Signal Processing, vol. 62, no. 12, pp. 3246-3260, 2014. 

[Sal88] D. J. Salmond, “Mixture reduction algorithms for uncertain tracking,” Tech. 

Rep. 88004, Royal Aerospace Establishment, Farnborough, UK, DTIC Number 
ADA197641, Tech. Rep., 1988. 

[Sal90] -, “Mixture reduction algorithms for target tracking in clutter,” in SPIE 

Signal and Data Processing of Small Targets, pp. 434-445, Los Angeles, CA, 
USA, 1990. 

[SKM95] D. Stoyan, D. Kendall, and J. Mecke, Stochastic Geometry and its Applications. 
John Wiley & Sons, Hoboken, NJ, USA, 1995. 

[SSS09] S. S. Stankovic, M. Stankovic, and D. Stipanovic, “Consensus based overlapping 
decentralized estimation with missing observations and communication faults,” 
Automatica, vol. 45, no. 6, pp. 1397-1406, 2009. 

[Str60] R. L. Stratonovich, “Conditional rnarkov processes,” Theory of Probability & Its 
Applications, vol. 5, no. 2, pp. 156-178, 1960. 

[SVV08] D. Schuhmacher, B.-T. Vo, and B.-N. Vo, “A consistent metric for performance 
evaluation of multi-object filters,” IEEE Trans, on Signal Processing, vol. 56, 
no. 8, pp. 3447-3457, 2008. 

[SWS09] D. Svensson, J. Wintenby, and L. Svensson, “Performance evaluation of MHT 
and GM-CPHD in a ground target tracking scenario,” in Proc. 12 th Int. Conf. 
on Information Fusion (FUSION), pp. 300-307, Seattle, WA, USA, 2009. 

[Tre04] H. L. V. Trees, Detection, estimation, and modulation theory. John Wiley & 
Sons, Hoboken, NJ, USA, 2004. 

[UCJ11] M. Uney, D. Clark, and S. Julier, “Information measures in distributed mul¬ 
titarget tracking,” in Proc. 14 th Int. Conf. on Information Fusion (FUSION), 
Chicago, IL, USA, 2011. 


140 



BIBLIOGRAPHY 


[UCJ13] 

[UJCR10] 

[VM06] 

[VSD05] 

[VV13] 

[VVC07] 

[VVC09] 

[VVP14] 

[VVPSIO] 

[Wil03] 

[WM06] 

[WvdMOl] 

[XBL05] 


—, “Distributed fusion of phd filters via exponential mixture densities,” IEEE 
Journal of Selected Topics in Signal Processing, vol. 7, no. 3, pp. 521-531, 2013. 

M. Uney, S. Julier, D. Clark, and B. Ristic, “Monte Carlo realisation of a dis¬ 
tributed multi-object fusion algorithm,” in Proc. 2010 Sensor Signal Processing 
for Defense Conference (SSPD 2010), London, UK, 2010. 

B.-N. Vo and W. Ma, “The Gaussian mixture probability hypothesis density 
filter,” IEEE Trans, on Signal Processing, vol. 54, no. 11, pp. 4091-4104, 2006. 

B.-N. Vo, S. Singh, and A. Doucet, “Sequential Monte Carlo methods for multi¬ 
target filtering with random finite sets,” IEEE Trans, on Aerospace and Elec¬ 
tronic Systems, vol. 41, no. 4, pp. 1224-1245, 2005. 

B.-T. Vo and B.-N. Vo, “Labeled random finite sets and multi-object conjugate 
priors,” IEEE Trans, on Signal Processing, vol. 61, no. 13, pp. 3460-3475, 2013. 

B.-T. Vo, B.-N. Vo, and A. Cantoni, “Analytic implementations of the cardi- 
nalized probability hypothesis density filter,” IEEE Trans, on Signal Processing, 
vol. 55, no. 7, pp. 3553-3567, 2007. 

—, “The cardinality balanced multi-target multi-Bernoulli filter and its im¬ 
plementations,” IEEE Trans, on Signal Processing, vol. 57, no. 2, pp. 409-423, 
2009. 

B.-N. Vo, B.-T. Vo, and D. Phung, “Labeled random finite sets and the Bayes 
multi-target tracking filter,” IEEE Trans, on Signal Processing, vol. 62, no. 24, 
pp. 6554-6567, 2014. 

B.-N. Vo, B.-T. Vo, N.-T. Pham, and D. Suter, “Joint detection and estimation 
of multiple objects from image observations,” IEEE Trans, on Signal Processing, 
vol. 58, no. 10, pp. 5129-5241, 2010. 

J. L. Williams, “Gaussian mixture reduction for tracking multiple maneuvering 
targets in clutter,” Master’s thesis, Air Force Institute of Technology, 
Wright-Patterson Air Force Base, OH, 2003. [Online]. Available: http: 
/ /handle .dt ic. mil/100.2/AD A415317 

J. L. Williams and P. S. Maybeck, “Cost function-based hypothesis control tech¬ 
niques for multiple hypothesis tracking,” Mathematical and Computer Modelling, 
vol. 43, pp. 976-989, 2006. 

E. A. Wan and R. van der Merwe, “The unscented Kalman filter,” in Kalman 
filtering and neural networks, S. Haykin, Ed. John Wiley & Sons, Hoboken, NJ, 
USA, 2001, pp. 221-280. 

L. Xiao, S. Boyd, and S. Lall, “A scheme for robust distributed sensor fusion 
based on average consensus,” in Proc. 4 th Int. Symposium on Information Pro¬ 
cessing in Sensor Networks (IPSN), pp. 63-70, Los Angeles, CA, USA, 2005. 


141 



